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SYMBOLS AND NOTATION 1 


A 


threshold friction speed parameter, u*.//p gD /p 
or the numerical factor in rarefied P ^ 

flows depending upon the nature of the reflection 
of the gas molecules from the particle surface 
(=3/1) 


A 1' A 2' A 3 ,A 4 


B 
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c 




D 

z 


e 



constant empirical numerical coefficients used in 
the friction threshold determination 

particle friction Reynolds number (=u*.D /v) 

L. p 

constant of the Riccati equation (=3K^ Pg/4D^p^) 

exponential constant used in correction factor of 
slip flow (=1.10) 

specific heat at constant pressure 
constant 

2 

drag coefficient [=drag force/ (1/2 pV S. ) ] 

J- a 

lift coefficient [=lift force/ (1/2 pV^S ) ] 

JL ri 

/ 

drag force 
vector drag force 
crater rim diameter 
particle diameter 

directional cosines of the vector drag force 

coefficient of restitution 

Eckert number (=U 0 /c (AT) 0 ) 

po 

function that relates vortex crossflow velocity 
to u*^ 

empirical wall effect function for vortex cross- 
flow 


Symbols and notation not in this list are defined and 
used locally within the text. 
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force 
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i'i/k 

k 

K 
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1 




Z 

£ t 

L 

L 

L* 


L x' L y' L * 

m 


m 


m 


apparent mass force 
Basset force 

resistance_ for. .molecular— regime- paxti'cl'e - motion ~ 

pressure gradient force 

acceleration of gravity 

mean velocity of gas molecules 

ripple height or reference height 

heat flux at planet surface or boundary layer 
shape factor (=6*/0) 

unit normal vectors of orthogonal coordinate 
system x, y, and z respectively 

von Kantian' s universal constant (%0.4) 

apparent mass coefficient, equal to 2/3 for solid 
spherical particles 

numerical constants associated with drag 
coefficients (see Appendix A) 

length 

wake length 

reference length or lift force ■ 
vector lift force 

3 

Monin-Obukhov Stability length, L*=u* c pT/kgH 
(where c is specific heat, T is ^ 

temperature, H is heat flux at surface) 
directional cosines of the vector lift force 
mass 


mass of gas molecules 
particle mass 



overturning moment exerted on a particle 
pressure 

Prandtl number (='c p 0 /y 0 ) 

Po 

magnitude of position vector or distance betwe< 
particle and vortex core 

position vector 

vortex core radius 

Reynolds number (^V^D^/v ) 

friction Reynolds number {=u*D^/v) 

crater rim height 

- 2 

bulk Richardson number (= (AT) 0 L 0 g 0 /T 0 U 0 ) 

2 

Reynolds number (=picD /v) 

P 

Rossby number (=U 0 /L 0 ft 0 ) 

Reynolds number (=V r D p /v ) 
shear (=3V r /3z) 

2 

surface area, equal to ttD /4 for a spherical 
particle ^ 

time 

temperature or time 

mean temperature at z = 0 

velocity 

final particle velocity 
gas speed 

Reynolds stress velocity terms 
freestream speed 
friction speed 


vii 



u 


*SL 


u 


*t 


U, 


*Un 


U,V,W 


U c 


u ,v ,w 

p p p 


V 


V 
-p 

V 


V 

-r 


V 


x 


V 


V 


V, 


V 


ef f 


w 

w 0 


W 


x,y,z 


x,y, z 


X 

^max 

z 

max 


local friction speed 
friction threshold speed 

friction speed of undisturbed velocity profile 

components of stream velocity x, y, and z, 
respectively. 

characteristic speed 

components of particle velocity x, y, and z, 
respectively 

stream velocity vector 

particle velocity vector 

magnitude of relative velocity ■ 

relative velocity vector 

x-component of relative velocity (=U-x) 

y-component of relative velocity (=V-y) 

z-component of relative velocity (=W-z) 

tangential velocity component in the vortex 

effect velocity component in the vortex — 

width 

\ 

initial or maximum vertical velocity of a 
particle's trajectory 

weight 

sample width 

coordinates of a particle's position 
components of a particle's velocity 
correction factor of rarefied flow 
maximum length of a particle's trajectory 
maximum height of a particle's trajectory 


viii 



Zo 

zj 

a 


6 


Y 

r 

AT 

(AT) 0 
6 

6 * 



e 


ijk 


n 


0 

e' 


K 

X 

u 

V 

p 

p. 

T 

<f> 


roughness height. 

roughness height in saltation 

fraction of gas molecules undergoing diffuse 
reflection at a particle's surface 

9 

slip coefficient (=r\/r\ ) 

0 

specific weight 

vortex circulation strength 

temperature difference. T-T^ 

temperature difference T 0 -T w 

boundary layer thickness 

boundary layer displacement thickness 

Kronecker delta (=0 for i^ j ; = 1 for i=j) 

alternating epsilon (permutation tensor) 

viscosity of medium 

coefficient of external friction 

number of gas molecules in one cubic centimeter 

potential temperature or momentum thickness 

fluctuation of temperature from the mean 

coefficient of thermal conductivity 

mean free path or ripple wave length 

absolute viscosity 

kinematic viscosity 

mass density of planet's atmosphere 
particle density 
shear stress or time 
dissipation function 


ix 



h'*2 

functions of the one-dimensional Riccati particle 
equation of .motion 

(0 

magnitude of rotation 

0) 

0, 

angular velocity vector . 

Q 

angular veloci^ 

Subscripts:^* 


( )o 

denotes' reference quantity 

< >m 

refers to model conditions 

i, j 

represents standard tensor notation 

A 

( ) 

nondimens ional quantity 

(") 

time averaged 

( )* 

instantaneous fluctuation from time average 

( ) 

vector quantity 


L Unless already presented 


the text. 


earlier or used locally within 
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I. INTRODUCTION 

The numerical calculation of particle trajectories under 
the influence of two- and three-dimensional turbulent boundary- 
layers in an incompressible Newtonian fluid flow and the 
experimental investigation of a three-dimensional downstream 
flow field around an idealized raised rim crater are the 
subjects of this study. Very few experimental investi- 
gations have been performed on three-dimensional turbulent 
boundary layer flow around a disturbing roughness element and 
none are known to the author with specific application to flow 
over a crater before this investigation was initiated. Also 
of interest was the numerical study of particle trajectories 
in a two-dimensional flow field under Martian surface condi- 
tions in which the effect of the rarefied atmosphere is 
included in the analysis. A combination of the particle 
trajectories and three-dimensional flow field around a crater 
is made by calculating a particle motion in the downstream 
flow field of a crater for Mars conditions. 

The study of particle movement and its associated 
trajectory is one with far reaching interests. An obvious 
case is that of wind erosion of soil and farm land. The study 
of snow storms, dust, and sandstorms and control is important 
to se’veral regions of the earth and an important problem of 
the future. Particle movement studies may even solve such 
distant and complex problems as the Martian eolian phenomena. 
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Also important is the understanding of the complex turbulent 
planetary boundary layer which moves the material. 

The trajectories and movement of particles by an eolian 
Jwind. blown) -process- ±s ~a~ complex "phenomenon. A classical work 
which deals with this problem is that of Bagnold (1941) in 
which he derives the relationship of the basic parameters 
based on data obtained from his field and wind tunnel studies. 
The yield of experiments Bagnold performed was equations for 
threshold friction speeds and mass transport of sand by wind. 
These relations were recently expanded from the specific case 
of movement of sand by wind to a general case of various 
density materials and any type fluid media motion by Iversen 
et al. (1973) and then to the case of low pressure flows (Mars) 
by Iversen et al. (1975a). 

Two-dimensional particle flow for surface conditions was 
calculated for both cases of ‘Earth and Mars. For the case of 
Earth, a turbulent boundary layer with a viscous laminar sub- 
layer and one without were calculated. For the case of Mars 
it was only necessary to calculate turbulent boundary layer 
flow with a laminar sublayer because of the low values of 
friction Reynolds number; however, it was necessary to include 
the effects of slip flow on a particle caused by the rarefied 
atmosphere. 

The effect of lift force on the initiation of particle 
motion has long been thought to- have a role in their movement. 
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This idea was analytically expressed by Saffman (1965, 1968). 

In the equations of motion, a lift force term was developed 
that acts on a single particle only in the laminar sublayer or 

i 

a corresponding small region of high shear near the surface for 
a fully turbulent boundary layer. The lift force functions 
were developed from the analytical work of Saffman (1965) for a 
single particle in simple shear flow. The lift force functions 
were consequently modified by empirical factors to account for 
wall effect and to match the limited experimental data 
available for Earth. An estimated interpolation for particle 
flow in the transition region for the case of Earth surface 
conditions was made from a combination of the laminar and 
turbulent sublayer solutions. The lift force, although rela- 
tively large near the surface, diminishes very rapidly, 
.increasing in height above the surface. 

With- use of the modified lift functions the simulated 
particle solution was numerically calculated for the Martian 
surface -conditions. A comparison is made between the effects 
of the surface conditions of Mars and that of Earth on the 
motion of the particle's trajectory. The effect of momentum 
loss due to an inelastic collision and rebounding of particles 
is considered for the Martian atmosphere. 

The advancement from a two-dimensional flow to the three- 
dimensional case significantly complicates the mathematical 
developments of the problem. The equations of motion of the 
particle become more difficult, but with the usage of high 
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speed computers can be solved numerically to sufficient degrees 
of accuracy. Once the equations are developed any type 
particle flow problem can be solved accurately if the flow 
f ield~i's" Known". Even with "just basic characteristics of the 
flow field known an approximate solution of particles can be 
obtained. 

The effect of turbulence (for Earth) on the particle's 
trajectories was modeled empirically by a cyclic fluctuating 
vertical velocity component and found to be of minor influence. 
The effect of turbulence on. the Mars trajectories was not 
accounted for after an analysis of its effect on Earth? cal- 
culation showed only a 4% or less altering of its trajectories 
for a particle of 100 microns diameter. Also the computing 
time to obtain a solution with the effect of turbulence 
included increased by an order of magnitude over the’ case with- 
out turbulence present. 

The path of the trajectory for a single particle was 
calculated for an idealized three-dimensional vortex system 
superimposed on a turbulent boundary layer flow under Martian 
surface conditions. Here a combination of the time— averaged 
velocity distribution and a Rankine vortex was made and 
particle .‘trajectories were numerically solved under its 
influence. An empirically developed function was used to 
simulate the wall effect on the Rankine vortex. The flow was 
assumed to be symmetric and stable around the model crater. 
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Several solutions of the particle's motion were calculated, 
each having a different relative position with respect to the 
vortex core. Two cases were calculated for the vortex motion, 
one of which was for a particle diameter of 100 microns and 
the other for a 500 micron particle. The results of the effect 
of a Rankine vortex present in a three-dimensional flow field 
simulating the wake crater flow show increased erosional 
properties outside of the horseshoe vortex system and also 
greatly increased particle trajectory heights. Also, smaller 
size particle trajectories followed the motion of the vortex 
flow causing lateral curvature of the particle's trajectory. 
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II. REVIEW OF LITERATURE 

What follows is a basic review of some of the aspects of 
the turbulent boundary layer. _ particle simulation-,- and- -crater 
modeling. It is not intended to be a complete or comprehensive 
compendium of these subjects. 

A. Turbulent Atmospheric Boundary Layer 

One of the basic features of a turbulent boundary layer 
is the turbulence. The most unique feature and the basic 
nature of turbulent motion is the fact flow parameters are not 
constant with respect to time at a fixed point in space, but 
fluctuate randomly through a wide .range of frequencies. The 
mean temperature and velocity profile are directly affected by 
the random fluctuation of temperature and velocity. Turbulent 
flow -has long been thought of as a three-dimensional flow with 
a random distribution of vortical eddies superimposed on the 
main flow. The mathematical expression for this was developed 
by Reynolds, i.e., 

u = u + u 1 (2.1) 

t 0 +T 

u » ~ | u dt (2.2) 

where u is the instantaneous velocity, u is the time averaged 
mean velocity, u' is the instantaneous fluctuating velocity. 
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t 0 is the initial time of the integral time period, and T is a 
sufficiently long length of time necessary to make u independ- 
ent of time. Of course there exist time averaging integrals 
for all variable flow quantities. By definition the time 
averaged integral of the fluctuating quantities is zero. 

The study of the turbulent shear layer has been wide. In 
particular the specific study of the turbulent boundary layer 
has been conclusive (Monin and Yaglom, 1965 and Lumley and 
Panofsky, 1964)-. From experiments conducted by Nikuradse 
(1932, 1933) on flow in rough-walled pipes, Schlichting (1968) 
and Sutton (1953) with others discuss the logarithmic wind 
profile law. For the case of a neutrally stratified planetary 
boundary layer Jensen (1958) and others observed that it also 
obeys the logarithmic law, i.e.. 


u(z) 



(2.3) 


where u(z) is the time averaged velocity which is a function 
of height, z; u* is the friction velocity? k is von Kantian' s 
universal constant; z is the height above the surface; and z 0 
denotes the equivalent surface roughness height. 

The two-dimensional turbulent boundary layer has been 
shown by Clauser (1956) to have a double structure. The 
double structured layer consists of an inner (or "surface") 
layer and an outer (or "defect") layer. The atmospheric 
boundary layer also has a similar type of double structure. 
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The surface layer is like the two-dimensional inner layer, 
while the outer layer is three-dimensional in nature. The 
three dimensionality of the outer layer is caused by a balance 

of- the— rotation - force "of" the 'earth (or planet) and frictional 

force. The resultant force caused by the rotation is called 
the Coriolis force. Coriolis force causes the direction of 
the mean wind in the outer portion of the planetary layer to 
turn to the right with increasing height (in the northern 
hemisphere) . The three-dimensional turbulent atmospheric 
boundary layer is also known as an "Ekman Spiral" or "Ekman 
layer" named after V. W. Ekman (1905) who first discovered the 
phenomenon and used it in a discussion of wind-generated 
ocean currents on a rotating earth. The rate of turning with 
height of the Ekman layer depends on the distribution of eddy 
viscosity and density. 


B. Laboratory Simulation 

An important problem in an experimental application of 
the atmospheric boundary layer is the laboratory simulation by 
use of wind tunnels. Great insight and understanding of the 
physical flow can be ascertained if correct similitude 
parameters are obeyed. A basic review of the criteria of 
correct model simulation is given here within. The review 
will be restricted to the lower portion of the turbulent 
boundary layer which, for the most part, is the most relevant 
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to the present problem. 

In order to properly model the lower portion of the 
planetary boundary layer in a wind tunnel certain necessary 
conditions must be met. 

Cermak et al. (1966), Hidy (1966), and McVehil 
et al. (1967) have shown that wind tunnels can model the 
atmospheric boundary layer. Many have applied the problem of 
modeling the atmospheric boundary layer by laboratory simula- 
tion with some degree of success (Halitsky, 1969; Plate and 
Ouraishi, 1965; and Cermak, 1963). Arya and Plate (1969) have 
shown that the Monin and Obukhov (1954) similarity theory 
offers a good foundation on which to base modeling of a stably 
stratified atmospheric boundary layer. , 

In the last fifteen years the study of surface (or inner) 
layer has been thoroughly extensive. As a result there are a 
large number of studies that have been conducted. Some good 
discussions of these are given in Monin and Yaglom (1965) , 
Lumley and Panofsky (1964) , and a basic review of the 
atmospheric layer in Monin (1970) . 

From the similarity theory of Monin and Obukhov (1954) 
the surface layer can be modeled if planar homogeneity exists. 
A detailed review of the implications of homogeneity to the 
similarity theory is given by Calder (1966). 

Exact modeling of the atmospheric layer in detail is not 
possible. However, by selecting certain similitude parameters 
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that need not be strictly matched and relaxing these require- 
ments, one can obtain good laboratory results for special type 
atmospheric flows. 

As reported by Cermak (1971) _-Uie_similari-ty- cri-feer-i-a—can 
be obtained from the equations of motion for the particular 
flow problem by nondimensionalizing the equations. As a 
result the nondimens ional equations will yield similitude 
governing parameters as coefficients of the equations of 
motion. 

i 

If horizontal and vertical geometry is kept it will result 
in an invariant nondimens ional transformation of the conserva- 
tion of mass'*' 


u. . = 0 and 


i£ + ll!V , 0 

9 1 9x . 


(2.4) 


Nondimensionalization of the time averaged momentum 

equation yields the criteria for dynamic similarity. The 

time averaged momentum equation is 

9u. _ 9u. _ 

9t + u j 9Xj + 2 e ijk fl j ^k 


1 __ 

P o 


8x i T 0 13 




9 2 u- ■ 9 (-u.u. ) 

i + 2 1 

3x k 9x k . 3x j 


(2.5) 


^Cartesian tensor notation and Einstein summation 
convention is used. 
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where the dependent variables are represented by mean (quantity 
represented by the bar over the variable and fluctuating 
values. The Boussinesq density approximation is made thus 
limiting the application of the equation to flows of AT << T„ 
where p is the deviation of pressure from the atmospheric 
pressure associated with p 0 . Using, 


A 

u. 

1 

= tu/Uo 

/V 

' u i 

= u|/u D } 

x i - V L 

r+ > 
II 

tUo/L 0 


= ; 

A. 

p* = p/pou: 


AT = AT/ (AT) „ ; g = g/g 0 (2.6) 


to nondimensionalize the momentum equation yields: 


9u. /s 9u. 

— + u. — r- 1 + 

3 


3t 


9x . 




u c 


A A 


2 e. ft. 


ijk H j *k 


(2.7) 


9x. 

i 


AT L ° g ° 

T„ U? 


AT g 6 i3 + 


u 0 L, 


a 2 ^ 

9 *k 9 *k 


9 (-u!u! ) 
9x . 


Three dimensionless parameters result as coefficients for the 
nondimens ionalization of the momentum equation. In order to 
maintain proper dynamic similarity the three similitude 
parameters must be matched. These three parameters are known 
as 
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Rossby number; R 0 =» u 0 /L o 0 0 

— - 2 

bulk Richardson number? IX = [ (AT) 0 /T 0 ] [L 0 g 0 /u 0 ] 
Reynolds- number'?' ' R~= "UpX o 7v 0 


The time averaged conservation of energy equation is. 


9T * n — 

at + i 3x i 


PoC 


2- 
3 T 




3 (-0 1 u! ) 

+ + 


PoC 


( 2 . 8 ) 


where <j> is the dissipation function. Nondimensionalizing as 
for the momentum equation it becomes 


3T . ~ 3T 

— + u i ~ 

3 1 3x d 


K 0 . 


V 0 

PoC p 0 V ° 


X* q XXq 

to J 


» — 


0 A 
3^T 




3 (-e*u£) 




r 2 

v 0 1 

Vo 

u 0 L 0 J 

c (AT) o 
_Po - 


+ etc. 


(2.9) 


Here two additional similarity parameters result: 

Prandtl number; P = v 0 p 0 c /k 0 

J- Po 

2 — 

Eckert number; E = u Q /c (AT) 0 

c Po 

if these two parameter criteria are satisfied then thermal 
similarity exists. For wind tunnel simulation in air the 
Prandtl number criteria is automatically satisfied and the 
Eckert number is only important in compressible flow. 
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For correct simulation of the pressure distribution (on 
various objects) in the atmospheric layer, in wind tunnel 
modeling Jensen (1958) observed that the roughness height ratio 
in the wind tunnel to that of the planetary boundary layer must 
be equal to a characteristic length of each, i.e., 




m 


( 2 . 10 ) 


Thus the geometric scaling of the boundary condition must also 
be satisfied as stated (Jensen, 1958) . If the other similitude 
parameters are also satisfied then dynamic similitude will 
prevail. A zero pressure gradient criterion in the longitudinal 
direction must also be satisfied. This can be accomplished by 
having a flexible ceiling of wind tunnel equipment and adjust- 
ing it until the criterion is met. In addition to the above 
criteria it is also necessary that a natural fully turbulent 
boundary layer be developed. This is usually accomplished by 
one of two methods, the first is a long test section in order 
to develop a boundary layer of the proper temperature distri- 
bution and turbulent structure. The second method is the use 
of turbulent tripping fencing and roughness elements. The 
latter is not as desirable for the reason that it usually does 
not obtain the characteristic of a fully turbulent flow 
because* basic flow variables change with distance downstream 
of the perturbating devices fairly rapidly. For the details 
of analysis see Cermak (1971) and Plate and Cermak (1963) and 



13 


Sundaram et al. (1972). 

In addition a few other basic points are pointed out. 

The geostrophic wind is not analogous to the free-stream 
.ve locity- -i-n-the-win d - tunn'eT. Wind "’tunnel studies are only 
valid for simulation of the surface layer of the planetary 
(or atmospheric) layer. Since this layer is essentially 
independent of the geostrophic wind and the Ekman spiral 
effect, it can be modelled without matching Rossby numbers. 

In order to simulate the entire Ekman layer the Rossby 
number criteria would have to be met which would involve a 
rotation of some nature of the wind tunnel test section. 

For the surface layer Monin and Obukhov (.1954) have 
developed a theory which uses only variable parameters defined 
at the ground surface, i.e., the surface roughness height, z 0 
and friction speed, u*, to entirely describe the layer. The 
surface layer extends up to a height of approximately 200 ft. 

Effects due to a nonneutral atmosphere (stable or un- 
stable) can be simulated to some extent by cooling and heating 
the wind tunnel floor and/or ceiling Cermak (1971) . 

One condition that must be met in order to obtain proper 
simulation is that the model flow be fully aero dynamically 
rough. This condition can be assured if the friction Reynolds 
number is greater than three, i.e., 


v 


> 3 


( 2 . 11 ) 
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For a detailed discussion of the limits of the wind • 
tunnel applications see Sundarum et al. (1972). 

1. Particle trajectories 

Very few scaled model experiments of sedimentation 
patterns around obstacles in air have been performed. The 
satisfying of the Froude number is not as crucial in air as 
compared to water, since in air there is no pertinent free 
surface. Of the experiments that have been , performed most 
deal with the accumulation of drifting snow. In one such 
study by Strom et al. (1962) the important similitude 
parameters were 

V Lo m' u 2 /g D p , V u ' e ' V u 

where L 0m is a characteristic model length, u is the stream 

speed at some reference height, D is the particle diameter, 

P 

Up is the particles terminal speed, e is the coefficient of 
restitution, and u^ is the speed of the snow particle. 

Other practical applications of the particle simulation 
study is that of the eolian phenomena that exists on Mars. 
Dust storms and other eolian activities have been suspected 
to occur on Mars oh the basis of telescopic observations 
(deVaucoulers, 1954; Kuiper, 1957? Rea, 1964; and others) and 
theoretical considerations of the Martian surface and 
atmosphere (Ryan, 1964;- Sagan and Pollack, 1969). Mariner 9 _ 
results confirm the existence of eolian features on Mars and . 
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show that eolian processes play a significant role in modifying 
the Martian surface (Sagan et al. , 1973). Many eolian features 
were observed in various stages of formation as the Martian 

dust storm- of _ i9 71-r9 72 ‘slowly subsided (Sagan et al., 1972). 

Analyses of Mariner 9^ imagery reveal many features that appear 
to have resulted from long-term eolian processes (Masursky, 
1973) . Knowledge of Martian eolian activity as a geologic 
process is essential to the understanding of the complex 
surface characteristics and geologic history of the planet. 

The modeling of the Martian eolian process uniquely lends 
itself to an investigation of the similitude parameters. 

The atmospheric pressure on Mars is much less than Earth, 
on the order of 5 millibars at the planet's surface. With 
these facts known, facilities with low pressure capabilities 
were conjectured to study the eolian phenomena (Bidwell, 1965 
and Chang et al. , 1968) . One such wind tunnel was constructed 
and low pressure tests were conducted to find threshold speeds 
at what was believed to be Martian surface pressure (Hertzler, 
1966a, b, and Adlon et al., 1969). For Martian threshold deter- 
mination it is necessary to use extremely low pressure, 
equivalent to the pressures on Mars, which can only be 
accomplished in a low density wind tunnel. However, the 
simulation of e'olian processes on the Martian surface can be 
accomplished by use of an atmospheric wind tunnel provided 
close attention is paid to satisfying appropriate modeling 
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parameters and proper simulation of the planetary boundary 
layer. 

Several important parameters of general particle salta- 
tion, sedimentation process have been found to be (Iversen 
et al., 1973): 

A Bagno Id* s coefficient 

B Particle friction Reynolds number 

Lift coefficient, C L = L/ (l/2pu^S A ) 

o 

C D Drag coefficient, C D = D/(l/2pu S^) 

D Drag force, (force) 

D c Crater diameter (length) 

0^ Particle diameter (length) 

e Coefficient of restitution 

2 

g Acceleration of gravity (length/time ) 
h Ripple height or reference height (length) 

i Length (length) 

L Reference length (length) 

Wake length (length) 

L* Monin-Obhokhov Stability length, L* = u* c pT/kgH 

ir 

(length) (where c is specific heat, T is temperature, 

P 

H is heat flux at surface) 

R Reynolds number, uL/v 

2 

S^ Reference area (length ) 

t Time ( time ) 

T Turbulence factor 

u Velocity (length/time) 
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u A Friction velocity, (= A/p (length/time) 

Threshold friction velocity • (length/time) 
u„ Terminal velocity (length/time) 

i? 

~a ~ Free stream velocity (length/time) 

Ug Geos trophic wind (above boundary layer) 

W Weight (force) 

W Sample width (length) 

s 

x Streamwise distance (length), 

y Lateral distance (length) 

z Vertical distance (length) 

z 0 Roughness length (length) 

z£ Roughness length of saltation (length) 

y Specific gravity 

6 Boundary layer thickness (length) 

k von Karman's constant (also C^) . 

X Ripple wave length (length) 

U Absolute viscosity (mass/length- time) 

2 

v Kinematic viscosity .(length /time) 

3 

p Mass density of atmosphere (mass/length .) 

3 

p Particle density (mass/length .) 

P 

2 

x Shear stress (force/length ) 

* These parameters can be arranged to form dimensionless 
similitude parameters which can be used to overcome some of 
the difficulties in modeling the full-scale conditions. 
Several important dimensionless parameters can be formed 
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e.g. (Iversen et al. f 1973), 

1. pD c /PpDp. By varying particle density and diameter and 
the crater diameter, this parameter can be varied from 
about 0.8 to 3. On Mars, this parameter would vary in 
value from about 1 for a 100 meter diameter crater to 100 
for a 10 kilometer crater. 

2. u{h)/u T _. Since the threshold friction speed u* is 

£ 

proportional to the reference velocity u(h) , providing 
geometry (including roughness) is exactly modelled, the 
ratio of reference velocity u(h) to terminal speed u^ 
will be modelled exactly if the ratio u*/u F is satisfied 
and if h/L is satisfied. 

3. [u(h)] /gL. The Froude number cannot always be satisfied 


in the wind tunnel without having a tunnel speed far 
below threshold speed. It is desirable to make' it as 
small as possible. Again since u* is proportional to 
u(h), this is equivalent to requiring a modelling material 
with as small a threshold speed as possible. The value of 
this parameter varies from 10 to 150 in the wind tunnel, 
and from approximately 20 for a 100 meter diameter crater 
to 0.2 for a 10 kilometer crater. 

e. The coefficient of restitution is* satisfied if model 


4 . 


and atmospheric materials have equivalent elastic 
properties . 
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5. S./L. Topographic features should be scaled exactly to 
satisfy this criterion. At large distances upstream from 
the region of interest, it is p rob ably_only_. necessary- to 
have equivalent scaled aerodynamic roughness. 

6. z 0 /L. The aerodynamic roughness should, in general, be 
to scale (Jensen, 1958) . Except for those craters 
surrounded by large-scale ejecta or other rough surface 
features ,' this is probably small on Mars. If the corre- 
sponding model surface in the wind tunnel is too smooth, 
it may be necessary to distort this parameter in order 
to obtain a turbulent boundary layer. It is important 

at the same time to insure that the ratio h/L be satisfied. 

7. zJ,/D c . If the equivalent roughness height in saltation 
z o is proportional to particle diameter, this parameter 
obviously cannot be satisfied on the laboratory scale 
model since such fine particles would have a very high 
threshold speed. Also, if introduced into the air stream, 
the particles would go into suspension and the saltation 
process would not occur. Calculations of saltation 
trajectory, however, show that the maximum height during 
saltation would be several times larger on Mars than on 
Earth, just as the saltation height on Earth is several 
times as large in air as it is in water. If the equi- 
valent roughness z' a is proportional to.uf/g, then z* 0 /D 

is proportional to Pp D p /pD c , the inverse of the first 
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parameter. 

8. h/L. The reference height h at which the reference speed 

is measured should be located within the logarithmic 

portions of the wind tunnel and atmospheric boundary 

layers. 

* 

9. z 0 /L . With a "naturally" developed boundary layer in 

the wind tunnel, a boundary layer velocity profile is 

achieved which corresponds to a neutrally stratified 

* 

atmosphere for which the Monin-Obhukov length L is 

« & 

infinite and the ratio z 0 /L is zero. A finite value of 
* 

L is achieved in the wind tunnel by heating or cooling 
the floor to obtain unstable or stable stratification. 
Another way of obtaining a nonneutral velocity profile in 
the wind tunnel (but perhaps not correct modelling of 
turbulence characteristics) would be by means of- shear 
fences, graded grids, or the like (Counihan, 1969). 

10. A/L. The relative ripple length may be related to zJ/L 
and the same comments apply. 

11. u„u*. and 12. u +J _D /v. As will be shown above, for a 

F *t *t p 

given condition such as for a modelling particle of 
diameter corresponding to minimum threshold speed, these 
two parameters would have the same values as for minimum 
threshold speed material on Mars. 

13. u*/u* t . The manner in which particles are transported 

and, in particular, the amount of material which is moved 
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is a function of this ratio. Thus, in order to keep u* 
as small as possible because of the Froude number, the 
threshold friction speed of the particle should be small. 

L-4-. u-(h') t/LT The Time^ c ale in the wind tunnel is much 

shorter than the time necessary for pattern development 
on Mars since the characteristic time is the ratio of- 
characteristic length L to reference velocity, u (h) . The 
time necessary for pattern development on Mars can thus 
be predicted from wind tunnel tests.- 

15. A Reynolds number u(h)L/v may or may not be an important 
modelling parameter.. For turbulent flows over sharp- 
edged features, the flow is - relatively independent of 
Reynolds number. The critical model Reynolds number 
(above which effects are independent of Reynolds number) 
depends upon model shape. If the model is too Stream- 
lined so that the test Reynolds number is below the 
critical, the model may have to be distorted by roughening 
the surface, creating sharper edges, etc. in order to 
lower the critical Reynolds number. Snyder (1972) quotes 
critical Reynolds number for sharp-edged cubes of 11,000 
and 79,000 for a hemisphere- cylinder . In the current, 
tests, Reynolds numbers based on crater diameter were 
generally above these values for sharp-rimmed model 


craters . 
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In addition to the previously mentioned criteria for 
similitude analysis further discussions can be found in Strom 
et al. (1962) and Warnock (1948). 

2. Crater modeling 

Martian eolian features occur in a variety of forms, most 
of which are associated with craters' or other topographic 
obstructions . The features are subdivided into two general 
types : dark streaks and light streaks (Figures 1 and 2). 

During the Mariner 9 _ mission, several areas were imaged 
repetitively in order to observe possible surface changes. 
Figure 1 shows a crater 17 km in diameter that .developed a 
dark fan-shaped streak within a 38 day period (Sagan et al., 
1972); it is typical of many dark features. Figure 2 shows 
light streaks associated with craters. These and all light 
streaks imaged repetitively during the mission showed no 
observable changes. Sagan et al. (1972) concluded that the 
light streaks are comparatively stable and the dark streaks 
unstable. Both types of features apparently ‘ can be used as 
surface wind direction indicators and some attempts have been 
made to derive surface wind patterns from streak orientation 
(Sagan et al. , 1973? Arvidson, 1974). 

A special interest case of this dissertation is the study 
of flow about impact craters under Martian surface .conditions, 
therefore a brief review is in order of the problems 




Figure 1. Crater 17 km in diameter in the Daedalia region near Solis Lacus ;on Mars, 
photographed on two revolutions (Rev. 1 15 on left. Rev. 195 on rilght) . 

The dark fan-shaped feature developed in a period of 38 days (from, Sagan 
et al. , 1972) • 1 


I 






Figure 2. Cratered terrain in Hesperia on Mars, photographed on two revolutions 
(Rev. 128 on left. Rev. 167 on right) , showing typical light streaks 
associated with craters and their apparent stability, compared to dark 
streaks (Figure 1) . Difference in contrast is the result of image 
processing (from Sagan et al . , 1972) 
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associated with wind tunnel simulation and also with the 
previous works already performed in the field. 

The physical simulation of the movement, depo sition and 

erdsTonTof f ine~ particles upon a complex terrain on a smaller 

geometric 'scale presents a complicated problem. The eroding 
of soil or sand is an intricate function of mean wind speed, 
frequency and intensity of wind gusts, particle size distribu- 
tion, density and shape of particles, surface drag forces, and 
the geometry of the topographic features. 

As Bagnold (1941) first reported there exists an optimum 
particle size that will correspond to a minimum threshold 
speed. The same universal trend is true for any planetary 
body. In the case of Mars, Iversen et al. (1973) and Greeley 
et al. (1973) report that both the optimum particle size and 
corresponding minimum threshold speeds are higher due to the 
low atmospheric density on Mars. The terminal speeds are 
approximately the same for like particles. 

The wind pattern over the crater would be similar to wind 
patterns observed over proturberances in boundary layers with 
small height to diameter ratios in laboratory scale tests on 
earth (Sedney, 1973; Gregory and Walker, 1951; and Hunt, 1971) . 
A horsehoe vortex is wrapped around the leading edge of the 
crater with the trailing vortices emanating downstream from the 
crater sides. The tangential component of velocity in each 
trailing vortex is outward away from the wake centerline near 
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the surface and inward above the vortex cores. The axial 
velocity components near the surface just behind the down- 
stream crater rim are minimal on the wake centerline with 
maximum shear stresses occurring on either side of the wake of 
greater magnitude than outside the wake {see Figure 3) . 

Further downstream the two surface shear maxima merge, and the 
maximum shear stress in the wake- is then on the centerline. 

Not all of the similitude parameters can be satisfied 
simultaneously in a scaled model experiment. In the case of 
modelling Martian craters the large geometric scaling factors 
which are necessary in order to properly simulate physical 
flow conditions cannot be met. The parameter D /L cannot be 

r 

satisfied in normal simulation facilities, e.g., if of a 
particle is 300 microns and- the crater to be modelled is one 

' _7 

kilometer, then the parameter D^/L would be 3 x 10 . This 

would mean that the particie size of a 30 cm model crater 
would be less than 0.1 of a micron. Obviously such small 
particles are not suitable since the cohesive and inter- 
particle forces would result in a relatively high threshold 
speed that when blown from the surface would go into suspension 
and therefore would not simulate the saltating phenomenon 
associated with the eolian process. A further complication 
that would result in the use of such small particles is that 
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FLOW OVER A RAISED RIM CRATER 



Figure 3. Raised- rim crater, showing horseshoe vortex. Axial 
velocity maxima (shown as vortex cores) of the 
trailing vortices converge downwind from the crater, 
forming a zone of higher surface stress than out- 
side the wake, and hence resulting in erosion 
(from Greeley et al. , 1974) 
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necessary, as in the case of sedimentation studies in water, 
to more effectively model the physical situation. 

The use of additional aides is necessary in order to 
determine the most important parameters or combinations of 
parameters of model simulation since all parameter similitudes 
cannot be established. One proven method used to gain addi- 
tional insight into the modelling problem is the nondimension- 
alization of the governing equations of the flow. Upon doing 
so several dimensionless parameters are manifested. After 
basic assumptions applied to the flow problem, the parameters 
needed to be' satisfied to ensure dynamic similarity are 
(Iversen, 1975b): 


V 

zj 


wind 

tunnel 



Mars 


( 2 . 12 ) 
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The satisfaction of Equation 2.13 is simple as it just 
requires geometric similarity of topographic and gross 
erosional and depositional features in the horizontal direc- 
tions . Small-scale- bed— form- features- such “as "ripple "wave- 

lengths would not be expected to scale without simultaneous 
satisfaction of all the original modeling parameters. 

The satisfaction of the other three similitude parameters 
is more difficult. Owen (1964) has shown that the equivalent' 
roughness height in saltation is approximately proportional to 
the ratio of the friction velocity squared to the gravity. 

Thus Equation 2.12 becomes 


Rh 2 gRH 2gRjj 
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*t 
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(2.16) 


but considering Bagnold's representation of the threshold 
friction velocity it becomes 
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Similarly for the Equation 3.15 
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(2.18) 


Assuming the threshold parameter A is a function only of the 
particle friction Reynolds number B, the drag coefficient is 
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(2.19) 


The parameter A is constant for larger particles and values of 
A and for earth and Mars would be approximately equal. 

This is true even accounting for recent evidence that B is not 
a unique function of A, even accounting for the rarefied 
atmospheric effects on Mars (Iversen et al., 1975b) . This is 
not true for small particles. The ratio of terminal to 
threshold friction speed u F /u* t is also a function of the 
particle and the range of values of this quantity on Mars can 
also be duplicated in the wind tunnel. The ratio of u*/u*^. 
is varied in the wind tunnel by adjusting the free-stream 
speed and so can be varied from values slightly below unity to 
several times unity as is also probably the case on Mars. 

Thus the groupings of parameters of most interest become 



pD c ®h 
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The approximate ranges of values of the individual parameters 
in the wind tunnel and on Mars are listed in Table 1. 

Table 1. Similitude__.par.ameter. .values 


Parameter 

Mars 

Wind 



tunnel 


Uj,/u* t 1 to 15 1 to 15 

u*/u*£ 0.8 to 2.0 0.8 to 2.0 

pD /p D 0.6 to 500 0.1 to 3.0 

cr K p p 

1/A 2 


50 to 100 


50 to 100 
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III. THRESHOLD DETERMINATION 

With the recent images received from Mariner 9^ showing 
evidence of eolian processes, new interest was generated in 
the conceptual understanding of the determination of static 
and dynamic threshold speeds under Martian surface conditions. 
This in turn created renewed interest in experimental and 
theoretical studies of the threshold phenomena under earth’s 
atmospheric conditions (Iversen et al. , 1973 and Greeley et al., 
1973} . On Mars the atmospheric pressure is approximately two 
orders of magnitude less than on earth and as a result the 
ratio of fluid to particle density on Mars is much less than 
earth's. This showed the need for experimental studies of 
particles with a greater density than a typical earth-like 
particle under earth’s atmospheric conditions and ultimately 
under Martian surface condition. 3- 

A completely analytical development of a threshold theory 
is very difficult and is usually amended to a semi-empirical 
theory with the aid of empirical constants found by experimen- 
tation. The following is a basic review of a recent paper by 
Iversen et al. (1975a) in which a combination of theoretical 
considerations and pertinent experimental data are combined to 

A low density atmospheric wind tunnel is presently under 
construction at NASA Ames Research Center, California. This 
facility would have the capability to duplicate the Martian 
atmospheric surface conditions. 
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develope threshold relations. The final form of threshold 
expression accounts for both Reynolds number and interparticle 
force effects. 

As reported by Bagnold (1941) there exists an optimum 
particle size which corresponds to a minimum threshold speed 
for a constant density material. The relation between the 
Reynolds friction number B and the threshold friction speed 
parameter is 


A = u* t / /p“gD^/p (3.1) 

B = u* fc D p /v (3.2) 

where u*. is the threshold friction speed, p and p are 
. P 

particle and air densities respectively, g is the gravitational 

I 

acceleration, D p is the particle diameter; and v is the 
kinematic viscosity of the fluid media. 

All variation of the parameter A occurs in flow where a 
laminar sublayer exists. A drastic increase occurs in the 
value of parameter A when the parameter B is less than one- and 
is decreasing. The experimental data of Bagnold (1941), ‘ 
Chepil (1945, 1959), Zingg' (1953), and Iversen et al. (1975a) 
support this claim. There is a large discrepancy in experi- 
mental data for friction Reynolds number less than five. 

Chepil (1945, 1959) explains the scatter by differences in 
particle size distribution. For small friction Reynolds 
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number B (very small particles) the parameter A is a function 
not only of parameter B but also of particle size and shape 
distributions as well as possibly other factors. Similar data 
scatter was experienced by White (1970) in threshold experi- • 
ments of small particles in liquid. 

Another factor is the cohesive properties of very small 
particles. The cohesive forces are caused by van der Waals* 
(London) forces, electrostatic charggs, and forces between 
absorbed film (Brown, 1961 and Kitchener, 1961). Particles of 
any size, if small enough, cohere on contact even when 
thoroughly dry and particularly in a vacuum (Gregg, 1961) . 

In metals, forces due to electrical interactions are more 
important (Kuhn, 1961) . The importance of cohesive effects 
for small particles was recognized by investigators interested 
in the erosive effects on lunar surface dust by lunar lander 
rocket engines (Roberts, 1966). Thus data scatter may be 
explained by the existence of interparticle forces due to 
moisture, electrostatic effects, and other forces of cohesion. 

The aerodynamic drag force and moment for a solid sphere 
resting on a plane surface and immersed in a uniform shear 
flow is calculated to be (Goldman et al., 1967 and O'Neill, 
1968) 

D = 1.7005 (6 7T p v S D 2 /4) (3.3) 

P 

M = 0.944 (2 it p V D 3 S/8) 

P 


(3.4) 
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Saffman (1965, 1968) derived an expression' for a transverse 
force in a simple shear flow for small friction Reynolds number 
(a first-order correction to the Stokes approximation) which 

.is, 

L = 6. 46p\>u(D^/4) / S/v (3.6) 

hr 

2 

For the case of the flow in the laminar sublayer u = u*z/v 
and assuming the center of the solid sphere to be at a 


distance of z = D /2 above the surface the transverse force 

P 

becomes 


L = 0.8077pu*D^(u*D /v) (3.7) 

Equating moments at threshold conditions about the down- 
stream point (s) of contact yields the friction threshold 
speed, however the moments due to drag, lift, weight, and 
interparticle force cannot be determined since they depend on 
the geometry of the sphere in relation to the surface of 
particles . 

At threshold the moments acting on the sphere are zero, 
thus equating these (lift, drag, weight, overturning moment, 
and interparticle force) and solving for the friction threshold 
speed results in 



1 + A 3 B 


(3.8) 
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where I represents the interparticle force and A^, and 

A^ are constant coefficients. The coefficients A^, A^, and 
A^ are unknown because the average geometry of the sphere's 
position relative to the surface of particles and fluid force 
coefficients are not generally known upon which these 
numerical coefficients are determined. An empirical set of 
numerical 1 coefficients are calculated by letting the 
coefficients A^, A 2 , and-A^ float for the data of Xversen 
et al. (1975a) and Weinberger and Adlon (1971) and solving 
for the least squares curve fit yields (Iversen et al., 1975a) 


u* t = 0. 266 


p gD. 
P P 


1 + 0.055/p gD^ 
P E. 


1 + 2.123B 


(3.9) 


This equation is the best curve fit to Iversen' s et al. (1975a) 

data points plus the four low air density data points of 

Weinberger and Adlon. The forms of the aerodynamic drag/ lift 

and moment are valid only for small Reynolds numbers (R = 

uD /v << 1 or B = 0.45). Equation 3.9 takes into account both 
P 

Reynolds number and interparticle force effects. The inter- 
particle forces are shown to be important in Iversen et al., 
1975a paper. 

Pollack et al. (1974) have shown that the corresponding 
threshold value of geos trophic wind will increase as non- 
erodible roughness increases up to a certain level and then 
will decrease as nonerodible roughness increases still further. 
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A very approximate expression has been derived which 
included both the effects of nonerodible roughness and inter- 
particle force. The expression for threshold friction speed 
is (I vers en-et-al; 7 — 1975 a 


a gD 

u * t ■ °- 257 /-V* 


1 + 0.0314 D °* 837 /p gD 2 
E P P 


l+0.221[iog e (l+D p /2z o ) 3 


(3.10) 


These equations should be regarded as tentative, however,, 
until further low-density experimental threshold data are 
available . 
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IV. PROBLEM FORMULATION' 

The type of flow where a particle is involved with a 
fluid media is called a two-phase flow system. When several 
particles are in motion the flow is termed multiphase flow, 
of which a two-phase flow system is a special case. There are 
basically three classifications of multiphase flow. The first 
is gas-particle flow, which has applications in fluidized 
beds, cosmic dust, nuclear fallout, and metallized propellant 
rocket flow. The second is gas-liquid drop flow with applica- 
tions in air pollution studies, gasoline and rocket engines, 
and many aspects of meteorological studies of the atmosphere. 
The third is liquid fluid-solid particle flow with applica- 
tions in studies of water sedimentation,, underwater machines 
analysis, and soil erosion by rivers and rainfall. The problem 
of interest here is solid particles saltating in air which is 
one specific case of the gas-solid particle two-phase flow 
system. 

As the flow of a gas over a surface of solid particles is 
increased from a very slow speed there occurs at a certain 
speed movement of particles which is caused by the net forces 
exerted on the particle by the fluid flow. These saltating 
particles are subject to three major forces, to the weight of 
the particle which tends to move the particles down toward the 
surface, a tangential force which maintains a forward motion, 
and. a normal force which is predominant near the surface and 
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moves the particle up from the surface and is caused by the 
pressure distribution on the individual particle's surface. 

The normal force is predominant near the surface due to the 
exis tence_o.f_.si:rong— shear —flow- between "the' particle In 
motion and the particles fo inning the surface. These forces 
exerted on the particle are resolved into two components, one 
parallel to the direction of mean flow, called the drag 
(tangential) force and the other normal to the flow, called 
the lift (normal) force (Raudkivi, 1967). Particle shape 
influences the magnitudes of these forces. The drag force 
consists of surface drag (viscous skin friction) and form drag 
due to the pressure difference in the front and back of the 
particle. The net drag force depends on the position of the 
particle when the force is applied and also on the position 
which the lift force acts when applied, both of these are 
functions of particle porosity, shape, size, and relative 
location of the particles (particle geometry) . In the follow- 
ing analysis the particles are assumed to be solid spherical 
particles of constant density. 

In addition to the three major forces there are several 
forces that are less dominant in this application. These 

l 

forces are the interparticle force, the force of the over- 
turning moment in shear flow, the Magnus force, the pressure 
gradient force, the Basset force, and the apparent mass force, 
and also the force resulting from the effect of temperature 
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gradients in the flow field. Not all the above forces apply 
directly to the saltation of particles being driven by a 
planetary boundary layer. It is necessary to examine the 
nature of each force and judge how much relevance it has to 
the specific problem of saltating particles. 

A. Drag Force and Coefficient 

Studies on the force that a fluid exerts upon a sphere in 
steady motion began with Newton's experiments in 1719. He 
found the force to be F, 

F = 0. 055irD^p (u-u (4.1) 

P P 

where (u-u ) is a relative velocity, D is the diameter of the 
p P 

spherical particles, and p is the density of the fluid medium. 
This relation seems to be fairly accurate for a Reynolds 
number region from 700 to 20,000 (Gugan et al., 1964), where 
inertial terms cannot be neglected as they are for Stokes flow 
in the momentum equation. 

Stokes reported in 1850 that for a symmetric flow field 
with small velocities (Reynolds number less than 1) , the 
dominant force on the sphere for these Reynolds numbers is 
viscous and the inertial terms in the momentum equation can 
thus be neglected. The Stokes force is (Lamb, 1932 and Stokes, 
1891) 
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F = 3irD y (u-u ) ('4. 2) 

P P 

where y is the kinematic viscosity of the fluid material. 

The drag fo 2 ^ce_is_ highly_ _dep.endent- -upon- the- Reynolds 
number. A standard way to manifest the Reynolds number 

i 

dependence on the drag force is to maintain the same form of 
the drag force equation for all range of Reynolds numbers. 

This can be accomplished by introducing a dimensionless param- 
eter called the drag coefficient C^, 


C 


D 


Drag force 
|p (u-u p ) 2 S A 


(4.3) 


where is a reference. In the case of a spherical particle* 
the appropriate reference area is S^, 

S a = ttdJ/4 (4.4) 

A p 

For Newton's experiments is a constant equal to 0.44. 
For Stokes flow C D is inversely proportional to the Reynolds 
number, the constant being 24, i.e., 

C D = 24/R (4.5) 

For flow where the Reynolds number exceeds one Stokes 
results are no longer valid. The reason is that the flow is 
no longer symmetric about the sphere because of the increased 
effect of the inertial forces. The resultant flow up to a 
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Reynolds number of ten has an oscillating effect which is 
difficult to express analytically. For flows' greater than ten 
boundary layer separation occurs - and a stationery ring vortex 
forms at the rear of the sphere. 

Oseen improved Stokes relation by including inertial 
terms in the flow field -behind the sphere. The drag coeffi- 
cient that results is a series expansion of the Reynolds 
number multiplied by Stokes relation; i.e., 

C D = iT [1 + yf- R + o (R z ) ] (4.6) 

This relation is approximately valid for Reynolds numbers 
up to approximately four. Proudman and Pearson (1957) 
improved Oseen 1 s relation by matching series expansion of the 
flow at the surface and away from the sphere. This resulted 
in additional terms of order Reynolds numbered squared 
multiplied by logarithm of Reynolds number (Prouridman arid • 
Pearson, 1957) . The region occupied by the vortex moves down- 
stream with an increase in Reynolds number until a Reynolds 
number of 150 upon which the oscillation of the vortex system 
begins . 

At Reynolds number roughly equal to 500 (in laminar flow) 
a wake pattern forms behind the sphere. Vortex rings are 
continually formed and shed from the sphere. This causes a 
periodic flow field which results in instantaneous unsteady 
values, of the drag force. 
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At Reynolds number greater than 30,000 the boundary layer 
on the spherical particle has changed from laminar to turbu- 
lent. The point of separation moves downstream due to the 
change from a laminar to turbulent boundary layer, this' causes 
a large decrease in the separated wake region behind the sphere 

which results in a significant drop of the pressure drag force 

\ 

of the sphere. This change accounts for a drastic reduction 
in the value of the drag coefficient. The value of Reynolds 
number at which the drastic drag decrease occurs is called the 
critical Reynolds number.' 

Free-stream turbulence affects the value of critical drag 
coefficient. Free-stream turbulence causing a. lowering of the 
value of the critical Reynolds number which causes backward 
shifting of the standard C D vs. Reynolds number curve. This 
results in effectively reducing the value of the drag coeffi- 
cient in the neighborhood of critical Reynolds number. The 
value of the critical Reynolds number was reported by Giedt 
(1951) to be dependent on turbulent intensity, characteristic 
length of the turbulent motion and the diameter of the spheri- 
cal particle. Generally, the drag coefficient is reduced due 
to the presence of turbulence and with increasing flow 
turbulence the critical Reynolds number continues to decrease. 

Since Newton's first experiments on the sphere, many 
investigators have experimentally developed the present C D vs. 
Reynolds number curve that is shown in Figure 4 (Schlichting, 
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Figure 4. The * standard' ' drag coefficient for solid spherical 
particles as a function of Reynolds number. Curve 
1 shows the Stokes theory solution, curve 2 
displays the solution to the theory of Oseen, and 
curve . 3 . records Newton*s estimate of drag 
coefficient (from Schlichting, 1968) 
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1968) . Also shown in the figure are the early drag coefficient 
relations of Newton, Stokes, and Oseen. This curve is valid 
for a single spherical particle moving in a continuum incompres' 
sible fluid at constant temperature and velocity. The effect 
of slip flow caused by rarefaction of the fluid will be covered 
later. The above description of flow about a. spherical 
particle for various regions of Reynolds number is covered in 
depth by Soo (1967) . 

In a- recent paper Bailey (1974) points out that many of 
the measurements of subsonic sphere drag have been affected by 
the experimental methods used in obtaining the value. His 
major point, among others, is that in the testing of spheres, 
if any support device is used it will usually result in a 
translational oscillation of the sphere which will increase 
the value of the drag. Thus he offers a revised "standard 
drag" curve to account for the inaccuracies manifested by the 
experiments. This revised "standard drag" curve does not seem 
to be needed in the calculation of saltating particle trajec- 
tories due to the fact that the Reynolds number rarely exceeds 
300 which is where Bailey's new curve begins to show variation 
from previous curves. 

For ease of computation of the drag coefficient calcula- 
tion, the method of Morsi and Alexander (1972) has been 
adapted and full details can be obtained in their paper. 
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The method Morsi and Alexander used is one where the 
curve of C D vs. Reynolds number is empirically fitted {to 
experimental points) in the form of 



K 1 K 2 
IT + ^ + K 3 


(4.7) 


The drag curve is broken into seven regions for the 
entire Reynolds number range. The length of each region is 
adjusted to’ negate any discrepancies at the end points. Stokes 
relation is used for Reynolds numbers less than one-tenth and 
a constant value of the drag coefficient of 0.4 is assumed for 
Reynolds numbers greater than 50,000. 

The method for obtaining the relationship of the above 
equation is the same for each region. It has the advantage 
over numerical techniques due to the fact that there exists 
an analytical solution for a simplified one-dimensional equa- 
tion of motion for a particle. If one-dimensional flow is 
assumed and the drag force is the only force on the particle 
the simplified equation of motion may be written as 


p dt 


- c d I o < u -y 2 s j 


(4.8) 


where m 
P 


is the mass of the particle given by 


m 

P 




(4.9) 


The above equation- may be rewritten in the form of a Riccati 
equation , i . e . , 
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du 

_J 

dt 


!>. - <j)„u + B,u 

1 2 p Ip 


where 


2 2 
3uk,u 3u k ? 3k,p.u 

h ” + ^' + "VT 

P P P P 


^2 = 


3yk 2 3k 3 pu 


. ^2 2p D 

4p D P P 

p p ^ F 


_ 3k 3 p 

B 1 = 4D"p” 
P P 


(4.10) 


(4.11) 


(4.12) 


(4.13) 


If u is constant an analytic solution exists, but if u 
is a variable the equation can be reduced to the Abel form 
that can only be solved analytically for special type velocity 
functions . 

Morsi and Alexander (1972) show the method for analytic 
solutions of various cases and go on to solve the constant 
coefficients k-^, k 2 , and k 3 for the seven regions of Reynolds 
number. In Appendix A the regions and values of the constant 
are given as found by Morsi and Alexander. The calculated 
values of C D are always within 2% of the experimentally found 
values of C D . 

With knowledge of the. characteristics of the drag 
coefficient the drag force term in the equation of motion may 
be developed. In a three-dimensional equation the drag force 
is 
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D = C r ,pV 2 D 2 TT/8 (4.14 

D r p 

where V r is the relative velocity vector of the particle 

V = V - V (4.-15 

-r - -p 

where V is the stream velocity vector and V is the particle- 
velocity vector, i.e.. 


V*» Ui+Vj+Wk 


where 


V = Ui + Vj+Wk 

-p p- p- 


U = x 
P 


V = y 
P 1 


W = z 
P 


V r = (U - x) i + (V - y) j + (W - z) k 


V r = / (U - x ) 2 + (V - y ) 2 + (W - z ) 2 


v r = |v r | 


where | V | is the magnitude of the relative velocity vector. 
The drag vector force may be considered as 


D = D i + D i+D k 
- x - v - z - 


(4.21) 
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where D , D , and D are the direction components of the drag 
x y z , / 

force. The drag force by definition always acts in the direc- 
tion of the relative velocity. ' Considering the situation in a 
two-dimensional case (omitting lateral motion of the particle) 
it can be easily seen that (see Figure 5) 


D =Ili D andD = ^D (4.22) 

x v r y v r 

where the magnitude of the normalized .relative velocity 
components are (U - x)/V r , (V - y)/V r , and in the three- 
dimensional case (W - z)/V r ; thus the x, y, z components of 
drag force for the equation of motion are: 

• 2 

x- component: V r (U ” X ) C D D p 71 ’/^ 

y-pomponent: V r (V - y)C D D^-n/8 

2 

z-component: V r (W - z)C D D^n/8 


B. Effect of Slip Flow on Drag 
Force and Coefficient 

The Knudsen number is the ratio of the mean free path 
between molecular collisions to a characteristic length of the 
flow problem (particle radius). In continuum flow the Knudsen 
number tends to zero but the Knudsen number increases as the 
flow becomes rarefied (hon - Maxwellian velocity distribution) . 
Slip flow results from the rarefaction of the flow and when 
the Knudsen number is 0.1 or greater slip flow must be 
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considered in order to keep the calculations valid. In the 
limit of Knudsen number tending to infinity the flow is called 
free molecular. The intermediate region is known as transi— 
tion flow for which very little theoretical work has been 

* »■, \ « t 

accomplished. The continuum and free ^molecular cases are well 
understood. 

In the derivation of Stokes relation it is always assumed 
that there is no velocity discontinuity at the sphere* s 
surface (continuum) . The relative velocities of. the fluid par- 
tide .next to the surface of the sphere are zero. .If there 
does exist a so-called “velocity slip" or velocity discontinuity 
at the surface of the sphere- then slip flow exists for which 
Stokes solution is no longer valid. An .approximate measure of 
the magnitude of the discontinuity is the product of the 
gradient of the velocity and the mean free path length of the 
gas molecule X. This discontinuity becomes noticeable either 
when X becomes comparable to a characteristic length of the ' 
flow problem or when the gradient of the velocity becomes 
unusually large (as is the case for a shock wave) . It is often 
assumed that the tangential forces acting on the surface of - 
• the sphere are proportional to the velocity jump, as described 
by Fuchs (1964), the resistance of the fluid to theparticle 
is given then by the relation 
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D = 37ryDp(u-Up) 


2 n + D p n e /2 

3n + Dpn e /2 


(4.23) 


where ri is the coefficient of external force factor of 
proportionality. Denoting the ratio n/h as 3/ commonly known 
as the slip coefficient the equation transforms into 


D = 6iry 


( u~ u ) D 

2 E 



1 +. 43/D 

P 

1 + 6 3/Dp 


(4.24) 


In the limiting case n e going to infinity the drag force 
returns to Stokes relation. 

Epstein in 1924 developed a rigorous analysis of the 
tangential velocity at the spheres surface and found 

3 = 0.7004 ^| - x] (4.25) 

where a represents the fraction of gas molecules undergoing 
diffuse reflection at the surface and (1 - a) the fraction 
undergoing specular reflection. 

For the case D p << X for very small particles or for low 
pressure conditions the motion of particles can be assumed to 
be of a molecular nature. If it is assumed that the particle 
creates no currents and the normal Maxwellian distribution of 
velocity is unaltered/ then the fluid resistance to the 
particle is a result of a greater momentum flux to the "front" 
side of the molecule than the "back" side. This is caused by 
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a greater number of molecules impinging on the front surface. 

r 

The resistance is proportional to the disturbance in the flow 
field which is proportional to the particle diameter squared. 
Then from statistical mechanics (for the mass of particle much 
greater than the mass of gas) the resistance is expressed as 


F = ? 6 n m G D iu - u ) 
m 3 g g g p P 


(4.26) 


where ri is the number of gas molecules/cm , G is the 
g g 

molecules mean velocity, (u - u ) is the relative velocity of 

P 

the particle, and 6 is the factor determined by specular 
reflection (6 = !)• or diffuse reflection (6 = 13/9) as recorded 
by Epstein (1924) . 

Again from statistical mechanics the mean free path is 
calculated by 


X = y/0.499n g m g G g 

then the fluid resistance is given by (D << X) 

(u - u ) 

F m = 2 + B) P_ 


(4.27) 


(4.28) 


Thus using the form of Equation 4.28 developed by 
Cunningham (1910) and applying Epstein's relation found in 
1924 yields. 
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6TT|iD (u “ U ) 

F _ p e! 

m 2(l+A2X/D p ) 


(4.29) 


where 


a = e/x 


(4.30) 


2 A 

Since 8/r = A — is small and A is approximately unity then 

p 

the above equation only holds for small A/D^. 

Knudsen and Weber (1911) performed experiments on 0.389 

cm radius glass spheres in reduced atmospheric pressure down 

2 

to 0.14 dynes/cm and developed a form of the slip correction 
equation for rarefied flow (noncontinuum) and subsequently, 
most of the "later investigators retained the form of their 
equation , i . e . , 


X=l + ~ [A + B exp (-CD /2 A)] (4.31! 

P - P 

where X is called the "Correction Factor" somewhat similar to 
the Cunningham factor. The correction factor divides the drag 
coefficient to yield the true result. Millikan in 1920 and 
1923 published the results of 10 years of experimental work 
and his values of A and B agreed within 1% of that by Knudsen 
and Weber. Millikan's experiment involved letting oil drops 
fall through air. Millikan's data should probably be 
considered the most reliable. Mattauch (1925) experimented 
with oil drops in nitrogen and had a region of Knudsen numbers 


from 0.1 to 5. 
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All of the above experiments were reviewed in 1945 by 
C. N. Davies who developed an average equation of all results. 
This lies very close to Millikan's results which are believed 
to be the most reliable. Figure 6 shows the comparison of the 
correction factor X of Knudsen and Weber , Millikan, Mattauch 
and C. N. Davies' averaged curve. The equation of C. N. 

Davies is used to calculate the slip flow effect along with 
the Chapman-Enskog calculation of mean free path. A compari- 
son was made and there was virtually no difference between 
Davies equation and Millikan's equations for X. C. N. Davies 
equation for X is 

,, -2.2D /X 

X = 1 + — [1.257 + 0. 4e P ] (4.32) 

P 

The above equation properly fits the two limits of 

X/D o and X/D -*• °°. 

P P 

In the hydrodynamic case of D >> X, the equation tends 

It 

to that of Cunningham's for flow very near the Stokes region; 
and in the molecular case X » D tends to the free molecule 

ir 

condition as described by Fuchs (1964). 

C. Lift Force and' Coefficient 

In view of many recent papers the lift force must be 
considered one of the pertinent dynamic forces that acts on 
small spheres in simple shear flow. For saltating particles 



A + B exp(-cD n /2A) 
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Figure 6. Slip correction for particle motion in gases 
(from Davies, 1945) 
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many investigators have shown that the existence of the lift 
force phenomenon (caused by the shearing flow) is necessary 
in order to initiate motion of the particles. 

Jeffrey in 1929 was one of the first to suggest that a 
fixed surface particle experienced a lift force caused by an 
unbalanced pressure distribution on the particles surface. 
Rubey in 1937 describes the force that moves a particle on a 
stream bed as the pressure difference between the bottom and 
top of the sphere among other methods. Then White in 1940 
reported that the lift force played a negligibly small role 
in the movement of particles by wind. This, the author 
believes was shown to be definitely false by several later 
papers . 

In one of these papers by Einstein and El-Samni (1949) 
the fluid dynamic forces on individual profusions from a 
hydraulically rough wall submerged in water were measured. 

The purpose of the paper was to study fluid dynamic 
forces acting upon particles in the surface layer of a sedi- 
ment bed over which a turbulent fluid flow existed. They 
measured the average 'lift force exerted away from the wall on 
0.225 ‘feet hemispheres that were glued to the wall in a 
hexagonal pattern. A lift coefficient could then be cal- 
culated as a result of the pressure differences that existed. 
The lift coefficient was found to have a constant value of 
0.178 if based on the velocity that was measured at a height 
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of 0.55 diameters of the sphere above the wall to which they 
were glued. 

Chepil (1958) performed a series of experiments similar 
to that of Einstein and El-Samni. It was found that the ratio 
of lift to drag force on the roughness elements at the surface 
remained nearly a constant value of 85% for a wide range of 
hemisphere diameters and friction velocities . A comparison 
is then made between Einstein and El-Samni 1 s results and • 
Chepil* s results under similar conditions. Chepil found an 
average C L value of 0.068. He then went on to point out that 
Einstein and' El-Samni* s results were based on a pressure dif- 
ference between the top and bottom of the hemisphere in 
contrast to his which used a more thorough pressure distribu- 
tion. Chepil also found that the pressure difference was 
2.85 times the lift that the hemisphere experienced, thus 
applying this correction to Einstein and El-Samni results 
yielded a C L of 0.0624. This value is only slightly higher 
than Chepil* s , results, who credits the discrepancy to the 
difference in the geometry of the roughness elements array. 

It may be of interest to note that in a recent paper by 
Willmarth and Enlow ( 1969) measurements were made on the 
fluctuating lift force acting on spheres at super-critical 

5 

Reynolds numbers (R - 4 X 10 ) and a fluctuating lift' 
coefficient of 0.06 was found. • The value of lift 

Jb 

coefficient is nearly equal to Chepil and the modified 
Einstein and El-Samni measurements, but Willmarth and Enlow *s 
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flow is laminar, while Chepil and Einstein and El-Samni's flow 
is fully turbulent, but both have relatively high Reynolds 
numbers . 

Chepil (19,59) in a later paper goes on to describe the 
effect of particle shape, particle density, angle of repose of 
the particle and closeness of packing of many particles as 
well as the impulses of fluid turbulence of the lift and drag 

I 

forces experienced by the individual particles. 

The discussion of the above material has been restricted 
to aerodynamically rough surfaces in contrast to a turbulent 
shear layer with a laminar highly viscous sublayer. The 
criteria or distinction between the two types of flow was 
first pointed out by Goldstein in 1938. He developed a fric- 
tion Reynolds number (formed in the usual way) based -on the' 
friction velocity of the shear layer, the particle radius, and 
the kinematic viscosity of the flow. He cited a critical 
value of 3.5 below which the laminar sublayer exists- and above 
which the flow is basically turbulent but in the process of 
transition. Schlichting later reported that this laminar- 
turbulent transition region exists to a friction Reynolds 

* t 

number of seventy, above which the flow is fully turbulent. 
Goldstein describes the viscous sublayer as a thin region 
having a high rate of shear which he believed was caused by 
the viscous stresses since the apparent Reynolds stress within 
the layer is small since the turbulent fluctuating velocity 
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components are negligible or nonexistent. Schlichting (1968) 
reports the layer as a laminar process where the viscous 
forces are much larger than the inertia force not allowing 
any turbulent fluctuations. 

Kline (1966) reports the viscous sublayer to be a time- 
dependent streaky structure that migrates slowly outward from 
the wall, and highly three-dimensional. Kline accounts for 
the streaky motion as a result of regions next to each other 
with relative velocities between them. A fairly regular wave 
pattern is formed by these, streaks in a transverse direction. 
Kline shows these transverse wave lengths to become shorter 
for adverse pressure gradients and longer for favorable 
pressure gradients. Kline et al. (1967) also observes that 
the unsteady three-dimensional streak flow of the sublayer as 
well as the turbulent boundary layer will disappear under a 
sufficiently favorable pressure gradient. This tends to 
support the claim that the viscous sublayer streaks are a 
phenomenon due only to the existence of turbulence above the 
sublayer. Kline's method of investigation was by use of 
visualization studies of dye and hydrogen bubbles of approxi-. 
mately 0.0005 inch in diameter which were supported by hot- 
wire anemometer data. 

Clark (1968) also conducted hot-wire anemometer tests 
that supported Kline * s results that the sublayer has an active 
role in the production of turbulence. 
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Bradshaw (1967) claims that local dissipation of . 
turbulence is. more than the local production but this is only 
true in the viscous sublayer and not true in the remaining 
inner region where they are equal. He then states that shear 
stress production is a result of active motion produced in the 
inner layer. Bull (1967) reports a similar type phenomena in 
support of Bradshaw's theory. 

Tritton (1967) takes an opposing stand that the streaks 
do not move outward thus causing a favorable effect for 
turbulent production. He deduces from measurements that the 
Reynolds stress makes a positive contribution at the wall and 
for streak flow to move away from or toward the wall would 
require a negative Reynolds stress contribution. According to 
Tritton in order to support the streak movement it would take 
a change in the signs of the correlation coef ficienti_of the 
axial and transverse velocity fluctuation which he states do 
not exist in his measurements. He explains that the streaks 
move away from the wall (in Kline's experiments) because the 
dye was injected in the stream thus giving them a velocity 
away from the wall. 

Thus, as shown by the opposing view presented, the 
viscous sublayer has some features that are still ambiguous 
and not uniformly accepted. 

The above discussion of viscous laminar sublayers has 
presented some of the basic features and nature of a single 
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phase (same fluid material) flow. With the induction of a 
single particle in the viscous sublayer the flow now becomes 
a two-phase system. It is assumed that the presence of the 
particle does not change the basic characteristic flow- field 
of the viscous sublayer. 

In a recent paper Bagnold (1973) reported several 
important results on the nature of saltation and bed— loading. 
One of these results was that Francis (1973) had demonstrated' 
that saltating of particles occurs in the absence of fluid 
turbulence in a laminar flow. This is an important result 
since in the absence of turbulence there exists no velocity 
components normal to the surface; thus if saltation occurs it 
must be predominantly due to the pressure distribution on the 
surface of the sphere, or a lifting force. Thus Bagnold and 
Francis have asserted the fact that the initiation of motion 
upward away from the surface is due to a lift force of the 
same order of magnitude as the drag force. Of course the 
existence of the lift force is short-lived as observation con- 
firms the rapid decrease of vertical force. 

Saffman (1965) derived an expression for lift on a sphere 
(computed to 0(R 1//2 )) moving through a simple shear flow and 
found 

L = 1. 615y (u - u p ) D 2 (S/v) 1/2 


(4.33) 
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where -S is the magnitude of the velocity gradient. 1 The above 
equation has the restriction: 


R << R 1//2 (4.34a) 

V K 

R « 1 (4.34b) 

1C 


where R v is the conventional Reynolds number based on sphere 

radius and relative fluid velocity (u - u ) . 

* P 

This relation can be utilized to analytically find an 
expression for the lift coefficient in a simple shear layer. 
For the case of the laminar sublayer the velocity profile may 
be written as 


u 


2 


u*z 


v 


(4.35) 


and the laminar sublayer height if z = 10v/u*, corresponding 
to a Reynolds friction number of 10 and gas velocity of u = 
10 u t . S is shearing rate. 


2 

3u _ U * 
9z v 


(4.36) 


1 In Saffman 1 s original paper the numerical coefficient 
was a factor of 4ir large. This error was pointed out to 
Saffman as a result of lift calculations and experiments 
performed by Harper and Chang (1968) ; Saffman (1968) then re- 
derived the equations and found an error in the coefficient 
of 4 it. 
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thus Saffman's lift formula for sphere in simple shear becomes 


L = 1.615y(u - u ) D 2 u*/v 

P P * 


(4.37) 


It is then assumed that the velocity of the particle in the 
initial stages of lift-off in the laminar sublayer is very 
small, and it is further assumed that the lift force only acts 
in the sublayer and possibly to some extent on the upward path 
of the particle trajectory, but is neglible in the downward 
part of the trajectory. Thus in the sublayer 


or the lift is 


u - u i u*z/v 


(4.38) 


L = 1.615pu*D 2 z/v 


(4.39) 


and the general lift coefficient for a spherical particle is 


p _ uij 

L irp(u - u p )^p 2 

therefore in a laminar sublayer C L becomes 

c = 12.923 v = 4.114 
U L it u*z .R^ 


(4.40) 


(4.41) 


Saffman's lift expression can be formed into a general lift 
coefficient expression for any shear flow as 


, = 4.114V 1/2 8u 

L 7u - u ) 3z 

P 


(4.42) 
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and is completely specified when the functional relations of u 

and (u - u ) with z are known, and should be approximately 
P 

valid if it obeys the small Reynolds number restriction. 

In the case of the aero dynamically rough surface (friction 
Reynolds number greater than 3.5), the above expression can be 
used with some validity. In the case of the rough surface 
there still exists a very high shearing region very close to 
the surface. This region of high shear rate is very thin and 
is similar to the viscous sublayer with the exception that the 
flow is turbulent. If the equation of mean velocity profile 
for rough surfaces is substituted in the above general 
expression for lift coefficient the result should be valid for 
an average lift coefficient as long as the restrictions 

R 2 « R « 1 (4.43) 

v K 

are not exceeded. Of course this represents a very thin 
region that exists from the surface roughness height z 0 at 
which the average velocity of the flow is zero to some finite 
height prescribed by the combination of the shearing rate at 
the surface, particle diameter, and kinematic viscosity of the 
fluid. The velocity profile for turbulent flow is given by 

u* 

u(z) = — log e (z/z 0 ) (4.44) 
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thus if the relative velocity is again to be assumed very 
small in comparison with the gas velocity 


u. 


y 


3u 


3z 


Kz 


(4.45) 


and the C L becomes 

C L = 3.186/[log e (z/z 0 )R f 1/2 ] (4.46) 

or letting z Q = D /30 
* P 

C L = 3.186/[log e (30z/D p )R f 1/2 ] (4,47) 

With knowledge of the lift coefficient expression the 
lift force terms of the equation of motion can be explicitly 
developed. The lift force is 

L = C T pirD 2 V 2 /8 (4.48) 

Li p J_ 

again where is the magnitude of the relative velocity 
vector (speed) . The lift- component of force acting upon the 
particle is in the positive normal direction of the relative 
velocity vector. The directional components of the lift 
force may be found by calculating the cross product of the gas 
velocity with the relative velocity and then taking that result 
and forming the cross product with the relative velocity again 
thus yielding the normal direction of the relative velocity 
vector in the direction of the lift vector. Then the direc- 
tional relationship of each component can be found by 
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normalizing the vector expression to unity: 

L | | V r x ( V r x V) (4.49) 

where L is the proper lift velocity direction but not the 
proper magnitude. Normalizing yields the proper direction and 
magnitude of the lift vector: 

L = [L i + L j+L k] L (4.50) 
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thus components of lift in the three-dimensional equation are 
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2 2 

x-component: L^C^pV^Dpir/S 

2 2 

y-component: L C_pV D ir/8 

y jj r p 

2 2 

z-component: L z C L pV r Dpir/8 

D. The Other Forces 

Several other forces act on the particle, few of which 
are of the same order of magnitude as the drag and lift forces? 
thus by recognizing certain basic underlying assumptions of 
the flow these minor forces can be shown to be of negligible 
contributions. One force, however, that is- of importance is 

, f 

the potential field that acts on the sphere and in this case 
is gravity. 

A force of considerable importance in calculation of the 
threshold speed at low atmospheric density is the interparticle 
force? but once the particle has initiated motion the force is 
of negligible effect in the calculation of the particle" 
trajectory. 

The force caused by the overturning moment is, as the 
interparticle force, important in the calculation of the 
friction threshold speed (Iversen et al., 1975a) as discussed 
in the previous section but not relevant to the particle 
trajectory motion, and thus is neglected in the calculations. 

The velocity gradient in shear flow causes the particle 
to rotate. For low Reynolds number this rotation of the' 
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particle causes a viscous interaction of the fluid surrounding 
the sphere. This is called fluid entrainment which effectively 
adds velocity to the one side of the sphere where the rotation 
direction is the same as the fluid velocity direction and re- 
tards the fluid velocity on the opposite side. This tends to 
move the particle in the direction of higher velocity, and is 
commonly known as the Magnus effect on spinning sphere 
(Goldstein, 1938) . Rubinow and Keller (1961) derived the 
following relation for the lifting Magnus force on a rotating 
spherical particle as 

l4d^p(mx (V - V )) [1 + 0(R)] (4.55) 

° ir ir 

where w- is the angular velocity vector of the rotating sphere? 
and the torgue on the sphere as 

L = -1 rrpD^ u[l + 0 (R) ] (4.56) 

r p 

which yields the equations of motion for the Magnus effect as 
I P D p -If- = - 3 ^ D p [1 + | R lY r + £ ^ p( s x V r > (4.57) 

1 dt = -ir ^ : p - (4.58) 

where I is the moment of inertia of the sphere. At high 
Reynolds numbers the separation point is shifted by the rota- 
tion of the sphere, this causes a force in the opposite 
direction as reported by Hoglund in 1962. This force results 
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for the particles whose diameter is smaller than the character- 
is tic length of turbulent eddies or the thickness of a shear 
layer (Soo and Tien, 1960) . 

The Magnus lift force is shown for a freely rotating 
sphere to = by Saffman to be an order of magnitude less 

than that caused by the lift force {resulting from shear flow 
without particle rotation) unless the rotation speed is very 
much greater than the rate of shear. In the case of saltating 
particles that are initially at- rest it is physically unlikely 
that they would reach rotation rates very much greater than 
the high shear rates of the lower most portion of a turbulent 
boundary layer whether or not a laminar sublayer exists. Thus 
in the trajectory calculations of spherical particle the ■ 
effect of the Magnus phenomenon has been neglected. 

Another term is the apparent mass associated with a 
sphere moving through a fluid. The fluid medium exerts a 
resistance against this co-called "apparent mass" as well as 
the sphere itself. The apparent mass of a body of revolution 
may be expressed as 

Hi = K' irR 3 p (4.59) 

ci 

where R is the radius of revolution and K* is a constant 
determined by the shape of the body, in the case of a solid 
sphere K '= 2/3. Thus the force for a mass constant sphere is 
given (of negligible magnitude in this case) 
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Next the pressure gradient force is considered. The 

pressure gradient 9p/9r is shown in Figure 7 for a particle of 

diameter D then 
P 




cos0 


lE 

[SrJ 


ttD 
■ P 


sin0 


D 

2 ^- d0 cos0 


(4.61) 


or 


F 

-pg 


„ irD t tt 

- iE —E f 

4 Jo 


cos 0 sin0 d0 


(4.62) 


The pressure gradient force after integration is 


~ 4ttD 
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(4.63) 


where p is the pressure (a scalar quantity) and r is the 
position vector of the particle. In the case of the planetary 
boundary layer it is assumed that planar homogeneity exists as 
well as a constant pressure with height through the inner 
layer of the boundary layer, therefore the pressure gradient 
force need not be considered. 

The Basset (1961) force accounts for the effects of the 
deviation in the flow pattern from steady state, or constitutes 
an instantaneous flow resistance of the particle. It was 
derived by Basset in 1892 as 
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where t 0 is the initial time when particle motion is considered 
and t is the final and instantaneous time. When the solid 
particle is accelerated at a high rate the Basset term becomes 
substantial as the observed drag force becomes several times 
the steady state drag force and the drag coefficient increases 
drastically. 

The force terms caused by the pressure gradient effect, 
apparent mass, and the Basset force are prevalent only if the 
density of the fluid is equivalent to or greater than the 
density of the particle. In normal application of saltating 
motion in the earth's and Martian atmospheres these effects 
are small in comparison to that of lift, drag, and grav ity and 
thus are neglected. 

The effect of the temperature gradient need not be con- 
sidered when the assumption of a constant temperature through- 
out the turbulent boundary layer has been made and it is 
assumed here. It is also assumed that no concentration 
gradients or radiation gradients exist. 

In Appendix B the resultant equations of motion are given 
for the two-dimensional and three-dimensional cases as well as 
the nondimensional equations for particle trajectory motion. 



74 


V. FLOW AROUND A CRATER 

The experimental determination of a three-dimensional 
turbulent boundary layer (shear) flow field in the presence 
of a perturbation element is one of the most formidable 
problems in current fluid mechanics research. The develop- 
ment of an analytical theory is also quite difficult since 
the effects of the strong viscous interaction in the neighbor- 
hood of the perturbation must be accounted for in three 
dimensions. Consequently there has been little progress in 
either the theoretical development or the experimental investi- 
gation. 


. A. Characteristics of the Flow Field 

The flow around an idealized impact crater may be 
considered to be a good example of the above situation. 
Basically there are two types of disturbances that exist in 
the flow field. One is with the presence of a so-called 
"large scale” disturbance element which is typically of a 
larger scale than the boundary layer height and the second is 
a "small scale” disturbance element. In the case of the 
"small scale ” , the disturbance of the flow field is contained 
within the boundary layer. The main difference between the 
large protuberance and a small one (relative to the boundary 
layer thickness) is that the small scale disturbance has only 
a local effect on the pressure gradient. In this case of flow 
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over a crater the disturbance may be considered small scale, 
and thus as Sedney (1973) describes the flow problem, in two 
parts, one in the immediate protuberance neighborhood and the 

other in the flow downstream. 

In order to have a good conceptual understanding of 
three-dimensionally disturbed flows there are several common 
characteristics to all flows whether or not the boundary 
layer is laminar or turbulent and regardless of the dimensions 
or geometry of the protuberance element and speed of the flow. 
In most cases the law of the wall for a boundary layer flow 
will break in the vicinity of the protuberance, but later 
downstream the laws will be valid when the effects of the 
disturbance have diminished. In the region of the disturbance 
the flow will experience streamwise vorticity (crossflow) . 
immediately upstream of the disturbing element one or more- 
vortices are induced. The primary vortex then stretches 
around the front of the disturbance (crater) and is termed a 
horseshoe-shaped vortex. The horseshoe vortex can be traced 
back to the secondary flow in the boundary layer upstream of 
the disturbance. A secondary set of vortices in an opposite 
sense of the primary pair is believed to exist on the outer 
side of the axial centerline next to the primary set. Another 
set of vorticies exists right behind the disturbing element 
that are closely spaced vortex filaments originating from 
spiral filaments which rise vertically behind the element. 
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The height of the filaments is approximately the same height 
as the disturbance element, but in contrast the horseshoe 
vortex is located closer to the surface. The sense of rota- 
tion of the horseshoe vortex is clockwise as looking down- 
stream from in front (upstream) of the element at the left 
hand side vortex which extends axially downstream. Gregory 
and Walker (1956) were the first to explore this type of 
phenomenon. These vorticies affect the velocity profiles of 
the flow by a redistribution of the momentum immediately down- 
stream of the element. The spanwise velocity profile in such 
cases was studied by both Tani et al. (1962) and Gregory and 
Walker. In the protuberance situation of three-dimensional 
flow in a boundary layer, vorticity stretching, concentration 
of vorticity upstream and downstream, and viscous effects must 
all be considered. Figure 3 displays a drawing of a horseshoe 
vortex system around an idealized crater model (Greeley, 1974) 

Tani (1968)- reports on examples of sudden perturbations 
given to a turbulent boundary layer such as roughness elements 
suction, or injections. He concludes that recovery to equi- 
librium is very rapid near the wall but rather slow in th'e 
outer region of the boundary layer. 

Although data are very limited on three-dimensional flow 
around a disturbance there are some experimental results 
available that support the existence of the. vortex systems. 
From the work of A. Hiderks, Prandtl (1952) presents pictures 
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of flow around a disturbance exhibiting the horseshoe vortex 
and two symmetrical spirals immediately downstream of the 
disturbance. Benson (1966) also pictorially shows the vortex 
system's existence for a hemispherical protuberance element. 

In conclusion to the brief introduction of flow about a 
disturbance it can be said that the resulting flow fields are 
complex but exhibit characteristic vortex patterns that are 
not unique to the flow conditions but are widely observed for 
many different flows. A very perceptive statement is found in 
Sedney's (1973) paper, "The three-dimensional perturbations 
found in experimental results are so complex in the neighbor- 
hood of the protuberance that it is unlikely that an analysis 
can be developed that is capable of describing that part of 
the flow field. There is hope that progress can be made 
towards analyzing the downstream flow field". 

This statement is found to be true by the present author 
from experimental investigations of flow around an idealized 
crater model. 


B. Experimental Results 

The main object of the experimental investigation of the 
three-dimensional flow field was to gain insight into the 
basic nature of the complex flow. A total description of the 
flow field was not sought but rather an in-depth study of the 
velocity distribution of the downstream flow for' the reasons 
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Figure 8. Wind tunnel model crater 
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mentioned above by Sedney (1973). with a complete knowledge 

of the downstream flow conditions a conceptual analysis of the 

intricate flow field is more readily obtainable. 

Two series of wind tunnel tests were performed in Iowa 

State University's wind tunnel. The tests were conducted in 

an open-circuit environmental wind tunnel which has an adjust- 

\ 

able ceiling in order to establish a zero pressure gradient. 
The test section is 6.5 meters long with a 1.5 square meter 
test section. 

The idealized model of the impact crater is shown in 
Figure 8. The model crater has a rim diameter of 30 cm and a 
surface to rim height of 4 cm. The height from the deepest 
point inside the crater (the center) to the rim is 8 cm. The 
overall diameter of the model crater configuration is 50 cm. 

The model crater was mounted on the axial centerline of 
the floor surface approximately 5 meters downstream from the 
entrance of the test section. This enabled the flow to 
develop as a naturally fully turbulent boundary layer in 
which the model was immersed. 

1. Pitot probe experiments 

The first series of experiments conducted used a vertical 
airf oiled-shaped probe. Pitot pfobes were secured in the 
center of the transverse width of the metal airfoil, thus 
aligning the axial direction of the probes with the freestream 
flow direction. This insured an exact measurement of the 
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axial component of the pressure distribution in the crater 
flow field. The probes were mounted on the airfoil to given 
equal values of the differences of the logarithm of the height 
above the surface for any two consecutive probes. This spaced 
technique is known as equal logarithmic spacing and is con- 
venient since the pressure and velocity profiles generally 
follow the logarithmic laws as earlier mentioned. This 
convenience manifests itself by displaying equal spacing in a 
semi -logarithmic plot of a height versus speed curve which 
will be used later. 

The vertical airfoil probe was then secured to a cross 
bar made of metal. This crossbar was mounted approximately 
one foot from the ceiling (so as not to disturb the flow on 
the floor below) on two metal grooved racks. The horizontal 
crossbar could then be moved in an axial direction. The 
vertical airfoil probes could be moved in the transverse 
direction. There were 24 pitot probes mounted in the vertical 
airfoil ranging from a minimum height above the floor's 
surface of 0.160 inches to a maximum height of 35.8 inches 
above the surface as shown in Table 2. This enabled the pitot 
probe system to measure the flow field three-dimensionally at 
any axial or transverse position and at 24 different height 
locations; the majority being close to the floor's surface 
because of the logarithmic spacing of the individual probes. 
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Table 2. Experimental pitot probe locations 


Position Height (in.) Log height 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 
16. 

17 

18 

19 

20 
21 
22 

23 

24 


0.160 

0.261 

0.332 

0.424 

0.541 

0.690 

0.881 

1.124 

1.434 

1.830 

2.335 

2.979 

3.801 

4.850 

5.635 

6.945 

8.550 

10.53 

12.915 

15.860 

19.440 

23.830 

29.210 

35.790 


-1.83243 
-1.58876 
-1.10140 
-0.85773 
-0.61405 
-0.37037 
-0.12670 
0.11698 
0.36067 
0.60433 
0.84801 
1.09169 
' 1.33536 
1.57904 
1,72900 
1.93802 
2.14593 
2,35423 
2.55839 
2.76380 
2.96733 
3.17095 
3.37451 
3.57767 
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The pitot probes were approximately 5 cm long with an 
outside diameter of 1 mm and an inside diameter of 0.75 mm. 

The pressure tubes that were attached to the. rear of the pitot 
tubes were polyurethane plastic tubing with an outside 
diameter of 2 mm and an inside diameter of 1 mm. The plastic 
tubing was carefully secured over the back end of the pitot 
tubes and then nearly laid against the side of the vertical 
airfoil section and secured then by fiberglass tape so not to 
disturb the flow field any more than was necessary. The 
plastic pressure tubes were run along the vertical probe to 
the horizontal crossbar and over to the metal grooved track 
along the side of the wind tunnel. They were then run through 
a hole in the surface of the adjustable wind tunnel ceiling. 
The hole was then patched with plastic wood to reduce any 
adverse effect on the wind tunnel flow. The tubes were then 
extended to a twenty-six port manometer board. Twenty- four 
ports were secured to the plastic pressure tubing and one port 
to a static port that was mounted in the wind tunnel in the 
vicinity in which the pressure measurements were to be made. 
The other port was left open to the atmosphere as a reference 
position. The maximum length of the movable fluid in the 
manometer board during the testing was approximately 60 cm. 

The manometer board tubes were made from 1/4 inch diameter 
glass tubes with an inside diameter of approximately 5/32 of 
an inch. The manometer board was inclined from the horizontal 
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upwards at an angle of five degrees. The fluid in the mano- 
meter Board was water of a soecial solution to retain good 
meniscus shape. 

All tests were conducted at an air speed of 75 ft/sec 
with an average deviation of ±2-3 ft/sec caused by different 
atmospheric conditions. The vertical probe was positioned in 
477 different locations where 24 different vertical height 
pressure readings were taken. A rectangular two-dimensional 
grid network was formed for one side of the axial centerline 
assuming the flow around the crater was symmetric. The 
measurements began at a position of ten inches axially down- 
stream on the centerline from the center of the crater model 
and extended downstream 12 inches giving 13 grid points with 
equal spacing of one inch. The measurements made in the 
transverse direction started at the axial centerline and 
extended outward 17.5 inches giving 36 grid points of equal 
spacing between points of 1/2 inch. In addition there were ■ 
nine separate measurements made on the transverse centerline 
starting at a distance of ten inches out from the center of 
the crater and extending to 13.5 inches out, again having 
equal spacing between grid network points of 1/2 inch. Thus 
the total number of two-dimensional probe locations was 477 
yielding 11,448 individual pressure points in the flow field. 
The majority of these measurements were made downstream of the 


crater model. 
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The, method of data reduction was to photograph the mano- 
meter board during the test and later to reduce the data from 
slides made of the photographs. This is believed to be the 
most valid method of. data reduction since it captures the. 
instantaneous location of all the meniscus positions and thus 
eliminates atmospheric and turbulent fluctuations that would 
occur if the measurements of all 24 probes were taken 
individually while the test was being conducted. Most of the 
fluctuation amplitude was damped by the long length of pres- 
sure probe lines. The manometer board was not photographed 
for each location until certain equilibrium was obtained among 
the meniscus tubes. The minimum time for this to occur was 
nearly five minutes but photographs were not taken until eight 
minutes to insure that equilibrium of the pressure tubes had 
been reached. 

The reduced data were normalized with respect to free- 
. stream to eliminate minor fluctuations caused by changing 
levels of turbulence in the atmosphere and minor variations of 
the freestream velocity. The data appear to be well behaved 
from this simple process. Thus all quantities (pressure, 
.velocity, etc.) were equal to unity at freestream conditions 
and the mean velocity is zero at the floor 1 s surface. 

The data obtained from the measurements were programmed 
on Iowa State University’s IBM 360/67 computer which has a 
simplotter graphics device capable of three-dimensional plots. 
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constant, line contour plots and two-dimensional plots. All 
477 two-dimensional surface grid network points were plotted 
(height vs. velocity) from the 23 individual height pressure 
measurements associated with each position. 

Figure 9 displays the coordinate system that will be used 
in reference to the model crater and later in the numerical 
study. The coordinate system used for the -model crater is 
right-handed orthogonal with the height above the floor's 
surface as the positive z direction. The positive y coordinate 
is in the downstream axial direction. The positive x 
coordinate is in the transverse (lateral) direction. The 
origin of the coordinate system is located on the vertical 
centerline of the crater at a height of 4 cm above the 
crater's floor at the mean level of the wind tunnel floor. 

Figure 10 shows a plot of the height z in inches above 
the surface vs. the dimensionless magnitude of the axial 
velocity component V (normalized by freestream speed) . Four 

V 

different experimental sets of data are displayed as well as 
the experimental curve of the undisturbed two-dimensional 
velocity profile for the turbulent boundary layer. The experi- 
mental undisturbed boundary layer thickness 6 was found to be 

I 

3.5 in. The model crater extended 45% into the undisturbed 
turbulent boundary layer. The four sets of data are 
individually taken along the constant position of y = 10 in. 
(immediately downstream of the crater) for x = 0,3,6, and 9- 
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Figure 9. Model crater coordinate system 
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Figure 10. Crater wake vertical velocity profiles for y = 10 in. 
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in. At the position x = 0 in. there is an unusually large 
velocity defect, almost one, very near the surface. There is 
a rapid increase in velocity with increase in height. At the 
lower end of the data points there appears to be a linear 
relation to log z with the dimensionless velocity and there 
also appears to be a similar relation for the large values of 
height and velocity. This is contrary to the single slope 
prediction of the undisturbed turbulent boundary layer. The 
slope of the linear line is equal to 0.4/u* from which the 
shear stress t can be deduced for the case of a zero pressure 
gradient. For the x = 0 curve this relation would not seem to 
be valid since it would have two slopes to choose from. This 
bi-linear relation also shows the possible nonequilibrium 
effect of the axial pressure gradient which would also violate 
the relation. 

As the x position is increased the velocity profile 
slowly approaches the shape of the undisturbed flow. As the 
downstream coordinate y is increased for similar plots the 
velocity defect becomes smaller and the undisturbed velocity 
profile is approached more rapidly with increasing transverse 
distance from the wake centerline. A very rapid acceleration 
of the flow' occurs over a relatively short distance. At a 
value of x = 12 in. only 2 inches further downstream, the 
value of the dimensionless velocity at z = 0.15 in. has 
increased from 0.024 to 0.52. As the distance downstream is 
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increased the velocity profile off of the wake centerline has 
a gradual transition from the characteristic turbulent boundary 
layer profile to a general type wake velocity profile and then 
returns €o the undisturbed profile far downstream. 

Figure 11 displays the height vs. dimensionless time for 
the case of a constant value of x = 0 in. for four different 
downstream positions, y — 10, li, 12, and 20 in. Along the 
axial centerline the axial component of velocity increases 
very rapidly with downstream distance near the surface. This 
is specifically highlighted by the drastic increase in 
velocity from the y = 11 to y = 12 in. position. Here the 
dimensionless axial component of velocity increases in the 
downstream axial direction. One interesting item is the shape 
of the velocity profile. As the distance downstream is 
increased the profile exhibits the characteristic shape of a 
laminar flow around a hemisphere with low turbulence levels 
with the exception of a velocity reinforcement downstream. 

The velocity reinforcement is a phenomenon wherein the wake 
speed near the surface is actually greater than the speed 
associated with the undisturbed flow at a similar height. As 
reported by Roberson and Chen (1970) and Rainbolt (1968) the 
strength and position of occurrence' of the velocity reinforce- 
ment is highly dependent on the level of turbulence in the 
flow. For high levels of turbulence the velocity reinforcement 
■if it occurs at all occurs sooner and is weaker than flow with 
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Figure 11. Crater wake vertical velocity profiles for x - 0 in. 
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less turbulence. 

Figures- 12-15 show plots of the dimensionless axial 
velocity component vs. the positive transverse length x at 
four different values of height z = 0.15/ 0.424, 0.881 and 
1.47 in. In each figure four different downstream values of 
y are presented, y = 10, 11, 12, and 20 in. as well as the 
corresponding value of dimensionless velocity associated with 
an, undisturbed turbulent boundary layer at the same height. 

All points in each figure at a constant height were measured 
by a single pitot tube and pressure tubing and thus any 
discrepancies that would be caused by use of different pitot 
probes are eliminated. Figure 12 shows the greatest relative 
velocity defect from that predicted by freestream. This seems 
plausible since it has the lowest value of height above the 
surface z = 0.15 in. A very rapid increase in the velocity 
defect occurs over a relatively short distance (2 in.) and 
then a slower gradual increase to stream conditions. Almost 
identical trends can be seen in the corresponding curves of 
Figures 12-15, with changes in amplitude. An apparently 
unique trend can be observed in all four figures at a position 
of y - 11 in. This trend shows a lower centerline velocity 
(x = 0 in.) and then a velocity spike at roughly x = 0.5 in. 
and then a return to the averaged increasing velocity defect 
curve. This velocity spike occurs at only the lower heights 
in the measurements and since the phenomenon was measured with 
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Figure 12. Crater wake lateral velocity profiles for z = 0.15 in. 
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Figure 13. Crater wake lateral velocity profiles for z = 0.424 in. 
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Figure 14. Crater wake lateral velocity profiles for z 
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different pitot probes it does not seem to be anomalous. If 
the flow is assumed to be symmetrical around the crater model 
then mirror images of velocity profiles would exist on the 
negative x side of the centerline. The resultant imagery of 
the velocity spike then matches the axial velocity profile of 
a horseshoe vortex system in laminar flow (Gregory and Walker, 
1951) . Although the flow in this case is turbulent a 
qualitative comparison can be made. The effect of turbulence 
on the horseshoe vortex makes it very difficult to obtain an 
instantaneous measurement of the flow with use of pitot tubes 
since they have the effect of time averaging the flow. In 
order to accurately observe such a vortex system the entire 
flow field must be instantaneously known, however such a task 
is not plausible from an experimental viewpoint. The meander- 
ing of the vortex cores makes measurements of the vortex 
properties impossible because of the time averaging effect. 

The geometry of the model crater to vortex system possibly 
makes the movement of the vortex small in the vicinity 
immediately behind the crater where the velocity spike is 
located. Although exact evidence of a vortex cannot be made 
from this analysis, the existence of some type of crossflow is 
shown. At a transverse length of 6 in. or greater away from 
the axial centerline the flow is within 90% of the’ undisturbed 
boundary layer velocity with the exception of Figure 12 where 
the velocity defect reaches out further transversely. As the 
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height is increased for the range of transverse direction to 
6 in. the disturbance of the flow downstream is greater than 
at lower heights (relative to undisturbed boundary layer 
velocity at that height) . This implies that turbulent mixing 
is greater downstream in the flow near the centerline. This 
region of greater mixing extends to a height of approximately 
the same as the disturbing element (crater). This doesn't 
seem unnatural as this also is characteristic of a wake. 

A constant velocity contour plot for a two-dimensional 
plane may be made with knowledge of the many velocity data 
points within the plane. Figure 16 is such a plot. The plot 
shows curves of constant dimensionless velocity (ratio of 
axial component to frees tream value of velocity) for the 
height vs. the transverse direction of the plane y = 10 in. 

A contour line represents an interpolated position for a 
constant value of velocity. In Figure 16 an increment of 
0.02 exists between contour lines.' The values of several 
contour lines are denoted in this figure and also in the 
following figures involving contour plots (thus enabling the 
value of the dimensionless velocity ratio to be known within 
two percent anywhere in the contour plane) . 

Upon examining the contour plot of Figure 16 it is 
evident that essentially for x greater than 8 in. the flow 
profile is logarithmic but has a much -slower speed near the 
surface than that of undisturbed corresponding flow (Figure 
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Figure 16. Crater wake lateral isotachs for y - 13 in. 
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17) for transverse distances greater than 8 in. Shear rates 
much greater than those associated with undisturbed flow exist 
for a range of height from 0.6 in. to 2.0 in. and extend 
transversely outward to approximately 6 in. (crater rim 
radius) . It may also be noted that there exist regions of 
very small velocity for a transverse range from the centerline 
(x= 0 in.) to x= 2.0 in. at heights less than 0.2 in. This 
seems to show the existence of the stagnant air region 
believed to exist immediately downstream of the crater 
(Greeley/ 1974). 

Figure 17 shows an identical contour plot of Figure 16 
with the exception that the plane is taken at an axial distance 
of y = 20 in downstream instead of y = 10. A comparison of 
these two plots shows much higher velocities near the surface 
and a drastically slower change of velocity with height at the 
x—z plane position of y = 20 than that of y = 10 in. A more 

regular pattern is formed in the y = 20 in. case that shows 

the flow is closer to equilibrium than the y = 10 in. case. 

Also in Figure 17 there are characteristics of wake flow near 
the axial centerline (x = 0 in. ) at heights less than 1 in. / 

e.g., the flow shows only a 6% increase on the axial center- 

line from a height of 0.1 in. to 0.6 in. 

At positive transverse lengths greater than 7 in. in the 
flow the velocity profile is very nearly that of an undisturbed 
flow. The velocity measurement of the pitot tube nearest the 
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Figure 17. Crater wake lateral isotachs for y = 20 in. 
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wall was never less than 68% of the freestream velocity or the 
flow again has obtained a high shear rate in the first 0.1 in. 
of height above the wall surface. 

The contour plots displayed in Figures 18 and 19 are of 
the y-z planes which are perpendicular to the x-z planes 
featured in Figures 16 and 17. Figure 18 exhibits the con- 
tours of constant velocity along the axial centerline (x = 

0 in.) from; y - 10 in. (located just behind the crater model) 
to y = 22 in. The plots show an extremely high acceleration 
of the flow along the axial centerline from y = 11 in. to 
y = 12 in. (also see Figure 10) . The speed of the flow 

measurement at the lowest height above the floor (z = 0.15 

/ 

in.) increases from 12% to 54% of freestream speed in 1 inch, 
then increases rather slowly with increasing distance down- 
stream. Again the flow at a height of 0.7 in. to 1.2 in. has 
a higher velocity gradient from y = 10 in. to y — 12 in. than 
the remaining downstream flow demonstrates. Figure 19 is 
similar to Figure 18 except that the position of the y-z plane 
is not at the axial centerline but instead at a position of 
x = 17 in. This contour shows very little influence on the 
flow caused by the presence of the crater. Sedney (1973) also 
noted that the disruption of flow caused by a "small" scale 
disturbance is most noticeable in the local neighborhood of 
the disturbing element and does not greatly effect the flow 
outside the immediate downstream vicinity of the roughness 
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Figure 18. Crater wake longitudinal isotachs for x = 0 in. 
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Figure 19. Crater wake longitudinal isotachs for x = 17 in. 
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element. 

The data that form two-dimensional contour plots may also 
be presented in a three-dimensional plot. Great insight of 
characteristic trends can be accompli shed^ by observing three- 
dimensional plots. Figure 20 shows a three-dimensional plot 
of dimensionless velocity vs. logarithmic height vs. axial 
distance. The logarithmic height scale is similar to that of 
the contour plots in Figures 16-19 . This is again used to take 
advantage of the logarithmic velocity profile of the turbulent 
boundary layer and also used to obtain relatively more data 
points in the lower heights where the velocity changes are 
greatest. Figure 20 is a three-dimensional display of the 
same data used to produce the contour plot of Figure 19. It 
again is along the axial centerline (x = 0 in. ) and extends 
downstream from y = 10 in. to y = 22 in. The coordinates of 
the four base points (A / B / C / D) of the three-dimensional dis- 
play are also given in the figure. Base point A (hidden in 
this particular picture) serves as the 'local' origin for the 
figure. The axial increments shown are in units of 1 inch and 
the height increments are in equal logarithmic height incre- 
ments. Several discrete values associated with grid network 
.points are displayed to give proper perspective to the figure. 
A close observation of this plot shows the same data displayed 
in Figure 19 only in a different fashion. The dramatic 
velocity increase can be easily observed in this plot. A 
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remarkable overall picture of the flowfield can be seen, with 
the net conclusion that the velocity defect caused by the 
crater is most apparent only in the first few inches immedi- 
ately downstream. 

A further comparison can be made between Figures 20 and 
21. Figure 21 displays the data three-dimensionally that was 
used to produce Figure 18. Again comparison between Figures 
18 and 21 will show ^jcact corresponding trends. The location 
of the y-z plane shown in Figure 21 is 17 inches off the 
axial centerline. Although this is only approximately 1.5 
crater diameters outside of the axial centerline, the flow 
remains nearly undisturbed except for a small velocity deffect 
immediately downstream of the crater. Again the y-z plane is 
located relative to the crater coordinates system by surface 
corner points A, B, C, and D. • • 

Figures 22 and 24 are corresponding three-dimensional 
plots of the contour plots in Figures 16 and 17, respectively. 
These plots are of the dimensionless velocity vs. the logarithm 
of the height vs. the transverse length x for x-z planes at 
constant values of y = 11 and 13 in. respectively. Figure 23 
shows an immediate position of y = 12 in. The large velocity 
defect caused by the presence of the crater can be observed. 
There is a rapid increase in velocity with height and 
transverse length. Again several discrete grid network values 
are displayed to reveal the perspective of the figures. The 
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Figure 22. Crater wake velocity surfaces for y = 11 in. 
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height increments are the same as in Figures 20 and 21 and the 
transverse increment is 1/2* in. An immediate velocity' defect 
can be observed in Figurp 23 and a very small velocity 
influence can be noted in Figure 24. By comparing the three 
Figures 22, 23 arid 24 which occur only over a 2 in. axial 
distance (y =• 10 in. to y = 12 in.) an overall detailed picture 
of the flpwfield changes both axially and transversely can be 
made. ' It might also be noted that three-dimensional plots of 

. y' 

similar nature to the above figures, only positioned further 
downstream, reveal no significant change over the figure 
presented in Figure 24. The overall velocity defect can be 
observed by a comparison of Figures 16-24. 

When the axial pressure and axial velocity profiles are 
known in the three-dimensional downstream flowfield of the 
crater a surface shear stress can be calculated from, empirical 
expressions . 

The empirical expression used to calculate the surface 
shear stress is one developed by Truckenbrodt (1955) for two- 
dimensional boundary layer flow at zero pressure gradient. 
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(5.1) 


where H is the shape factor defined as 6 /0,5 is the boundary 
layer displacement thickness and 0 is the boundary layer 
momentum thickness defined in the usual manner. 
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Although the first couple of inches of the measured 
velocity field may not be in a zero pressure gradient flow, 
and the flow is not two-dimensional, some insight can be gained 
about the surface shear stress distribution by examining cal- 
culated stresses and the trends they exhibit. Each of the 477 
velocity profiles measured were numerically integrated to 
calculate, the local displacement and local momentum thick- 
nesses. With knowledge of the individual values of displace- 
ment and momentum thicknesses the surface shear stress can be 
calculated from Truckenbrodt ' s empirical relation from which 
the local value of the friction velocity u*^ may be determined. 
Figure 25 shows a contour plot of the constant ratio of local 
friction speed to the friction speed associated with un- 
disturbed flow under identical conditions for the axial 
distance downstream y vs. the transverse length x. The 
constant contour lines are in increments of 0.10 of the 
dimensionless friction speed ratio u^/u*. 

Upon examination the contour lines of Figure 25 show a 
region of high surface shear stress (friction velocity) that 
lies immediately downstream and in the wake of the crater. 

The region of high friction velocity (predicted by Equation 
5.1) unfortunately may be invalid but the trends displayed by 
it seem to be correct. A comparison may be made with experi- 
mental data taken by Udovich (1973) . Udovich performed a 
series of tests on the turbulent boundary layer flow around a 
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singular idealized crater model. The crater had a 5 in. rim 
diameter and a 7 in. overall diameter. The tests conducted 
were to place a 1.5 cm circular patch of erodible material at 
only one location in the crater flow field per test. The 
freestream was then increased from a low value to a freestream 
velocity that was sufficient to move the material within the 
circular patch. Forty-five various positions were measured 
twice and recorded. It is assumed that the local freestream 
speed is proportional to the local friction speed. Then the 
ratio of freestream speed necessary to move material in the 
absence of the crater ,to that necessary to move material in 
the presence of the crater is equal to the ratio of local 
shear stress to the local shear stress in the absence of the 
crater, i.e. 
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Figure 26 displays the results of Udovich's data if the flow 
field is assumed to be symmetric and the mirror image posi- 
tions on opposite sides of the axial centerline are averaged. 
It is noted that there exists a large amount of scatter for 
each data point. Figure 26 uses the same coordinate system as 
Figure 25 with the exception that it is nondimensionalized by 
the crater rim diameter. The dashed line represents the' 
experimental results determined on the curved slope of the 
raised rim diameter. The shaded area of Figure 26 represents 
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the same area presented in Figure 25 if each are compared on 
a dimensionless basis. A comparison between the shaded sec- 
tion of Figure 26 and Figure 25 shows similar qualitative 
trends for the ratio of friction speeds but quite a difference 
in values of magnitude. For the regions of high surface shear 
stress the true value of the friction speed ratio should be 
somewhere between the two values given in the two figures, 
since Udovich's experimental technique has a tendency to 
underestimate the true value of the local friction speed and 
Truckenbrodt' s empirical results tend to slightly overestimate 
the local friction speed. The numerical results obtained from 
Truckenbrodt * s relation should be assumed to be more accurate 
than the data of Udovich. This may possibly explain the large 
differences in the ratios of u*^/u A of Figures 25 and 26. 

The value of the friction speed ratio may also be a function 
of the Reynolds number if it is below the critical value. For 
Udovich's experiments the value of Reynolds number is just 
below the critical value; however, the value of Reynolds 
number for the experiments performed in Figure 25 was above 
the critical value. This too could possibly account for some 
of the differences. 

2. Hot-film experiments 

The second series of experiments conducted on flow around 
the model crater were made with the aid of a hot film velocity 
measuring device. A Thermo-Systems Inc. model 1080 total 
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velocity vector anemometer and model 1126 calibration system 
were use^i to measure the three-dimensional velocity vector. 

The system 1080 measures the velocity magnitude and direction 
over the full 360° solid angle flow fields. The device is 
capable of measuring speeds up to 300 ft/sec at one atmosphere 
of air and claims an accuracy of ± 3% of the magnitude of each 
direction ± 0.1% in total magnitude. It further claims two 
independent directional angle measurements each with 3° over 
the complete solid angle (4 tt Steradians) . The temperature 
range is from 0°F to 200°F with capabilities of measuring the 
temperature within ± 2°F. The standard time constant of the 
electric circuitry is 20 milliseconds at 60 ft/sec at standard 
conditions although this can be altered by use of an active 
filter. The frequency response is one kilohertz direct 
current . 

The system includes a streamlined probe which houses the 
orthogonal triad of hot film sensors, a 15 foot power cable, a- 
power control circuit, a supplementary active filter, and also 
the pneumatic calibrating system. The three orthogonal sensors 
are single-ended rods which have a diameter of 0.006 in (0.15 
mm) and length of 0.08 in (.2 mm). - Each sensor rod is composed 
of a quartz rod with a thin conducting film made of platinum 
0.0002 in thickness covered by a patented quartz coating. The 
data output of each measurement consists of six analog voltage 
signals, each of which has a possible range of 0-20 volts. The 
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output can either be measured by a digital voltmeter or dis- 
played on an oscilloscope. The data can also be stored on a 
tape recorder or oscillograph. 

The method used in these experiments was to use an active 
filter with a time constant setting of 30 seconds for each of 
the six voltage signals. It was necessary to time average 
each signal since the entire flow field could not be measured 
simultaneously and in order to compare results an average 

was needed. Time averaging has the tendency to wipe out 
details of the vortex systems caused by meandering of the 
vorticies. However, large scale trends relative to the 
boundary layer thickness should still be observed. 

The model 1080 probe was mounted on the vertical airfoil 
probe in similar fashion to the earlier described test. The 
cable leading from the rear of the probe was carefully secured 
above the surface downstream to the side of the wall and 
through the corner of the wind tunnel wall and floor to the 
control panel. 

Experiments were made at two different transverse planes 
with four height measurements in each plane, 3/4, 1-1/2, 

2-1/2, and 3 inches. These measurements were made in order to 
determine the lateral, vertical and axial speed components. 

The axial lengths downstream were y = 12 in. and y = 18 inches 

f 

for the two transverse planes which were measured. Thirteen 
transverse positions were measured at each height and each 
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plane yielding a total of 104 individual positions of three- 
dimensional velocity measurements. Each measurement took 
approximately 45 minutes to make.- Before each session of 
testing the model 1080 system was calibrated and adjusted if 
necessary to identical settings for all measurements made. 

The frees tream speed was approximately 50 ft/sec ± 2 fps/sec 
in order to obtain greatest accuracy as described in the 
Thermo-Systems Inc. Instruction manual (1970). The data were 
recorded and reduced according to data reduction procedure as 
.found in the Instruction manual. A data reduction technique 
described by Teilman et al. (1971) as a so-called 'improved' 
technique was also used for various positions and no signifi- 
cant difference existed. 

Figure 27 shows a plot of height vs. transverse length 
for an axial location of y = 12 in. adjusted to s ymme tric flow. 
At each position where a measurement was made a dimensionless 
velocity vector is plotted for the crossflow components w and 
u. The velocity vector is made dimensionless by dividing the 
crossflow magnitude by friction speed associated with free- 
stream flow. Several values of the crossflow velocity ratio 
are given to show relative magnitudes. Within a transverse 
length of 6 inches there appears to be a general turning of 
the crossflow inward and also upward. This appears to be a 
plausible result since the location of the z-x plane is at 
y *s 12 in. whereas the flow was earlier shown to be in a high 
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Figure 27. Crater wake crossflow velocity vectors for y - 12 in 
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velocity gradient region. The flow must return inward down- 
stream in order to form the high velocity of the wake region 
downstream of the low-speed region formed immediately down- 
stream of the crater. The absolute time-averaged magnitude of 
the crossflow at any position measured did not exceed 7% of 
the total magnitude of the velocity vector. The magnitude 
(strength) of the crossflow would of course depend upon the 
geometry of the disturbing element that gives rise to the 
formulation of the vortex (crossflow) system. The smaller 
scale effects of the vortex system would not be expected tp 
be observed in this type of time averaged measurement due to 
meandering of the vortex cores . 
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VI. NUMERICAL SOLUTIONS 

The equations of motion for a spherical particle under 
various given conditions were solved by a numerical computing 
technique which employs the use of a scientific subroutine 
package called node. Node calls six other subroutines, three 
of which must be supplied by the user. The solution of first 
order ordinary differential equations with initial conditions 
is obtained by using the predictor-corrector equations of 
R. L. Crane (1962) . These equations have a wide range of 
stability. The necessary back points are initially calculated 
with the Ruhge-Kutta-Gill single-step method. The corrector 
procedure is not iterated. Node has the ability to automatical- 
ly check the solution's accuracy at each step and change the 
stepsize to meet the specified accuracy of the user if it 
exceeds the error limit. It also will double the stepsize 
(provided accuracy is met) in order to reduce computing time. 
Copies of the computer program used for the solution of the 
particle trajectory equations of motion may be obtained at the 
Department of Aerospace Engineering, Iowa State University, 

Ames, Iowa, 50010, with reference made to this study. 

All particle calculations have several common elements 
involved in the numerical solution. Table 3 displays the 
values of the constants of the acceleration of gravity , 
coefficient of the kinematic viscosity, and .the density for 
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Table 3. Calculation 

variables 


Variable 

Earth 

Mars 

2 

g (ft/sec ) 

32.1725 

12.5679 

v (ft 3 /sec 3 ) 

1.5 x 10“ 4 

116.786 x 10 -4 

p (slugs/ft 3 ) 

2.25 x 10" 3 

2.00 x 10~ 5 

p (slugs/ft 3 ) 

br 

3, 5, and 7 

3, 5, and 7 

u # (fps ) 

1,2,3, and 4 

10, 15, 20, and 25 

0^ (microns ) 

50 - 1050 

50 - 1050 


both Earth and Mars. Also shown in the table are the various 

ranges of particle density p , friction velocity u*, and 

P 

particle diameter D for which particle motion was calculated. 

P 

The numerical computing scheme specified four place accuracy. 

Most particle trajectory calculations were initiated with 
a value of height one-half diameter above a plane surface. 
Since a particle's geometric relation to an actual surface 
composed of like particles is not unique but rather is one of 
a random distribution of heights, the calculation is assumed 
to represent an average particle trajectory. Several 
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different initial values of height were calculated with many 
different types of lifting functions (all functions of height 
z) . These calculations exhibit a very nearly constant 
particle motion for initial values of height of three-tenths 
to seven-tenths of the particle diameter. Above and below 
these values, particle motion was drastically curtailed. A 
possible explanation of this is that in the case of initial 
values of height less than 0.3 the air speed is much less than 
that of the higher heights and thus the lift force is much 
smaller. In the case of initial values greater than 0.7, 
since the lift coefficient C L is a function of height and 
becomes very small at larger heights, again the lift force is 
small. 

All values of drag coefficient C D are calculated from 
the equation listed in Appendix A. For the case of the 
rarefied atmosphere on Mars a correction factor based on C. N. 
Davies (1945) equation is applied to the equations of motion 
for a value of mean free path, X = 13 microns. The mean free 
path of the Martian atmosphere was calculated from Equation 
4.27 assuming an atmospheric surface pressure of 5 mm of 
mercury. All particles were assumed to be point masses 
located at the center of the particle. 

The coordinate system used for two-dimensional particle 
trajectories is one for which the positive z-direction, 
measured from the mean level of the surface, is upward away 
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from the surface. The x- direction is aligned in the downstream 
direction of the boundary layer flow. For the case of three- 
dimensional simulation of flow around a crater the coordinate 
system is that described in Section V with' corresponding 
Figure 9 . 

Figure 28 shows a plot of the particle friction Reynolds 
number versus size of the particle diameter for Earth and Mars 
cases. The solid lines denote the Earth case for the dif- 
ferent values of friction speed u* and the dashed lines 

CO 

represent Mars. As can be observed from this figure all cases 

of the calculations for Mars lie in the regime of a turbulent 

boundary layer with the existence of a laminar sublayer. The 

majority of the Earth cases calculated lie in the transition 

region with the exception of the smaller diameters for each 

friction speed for which a laminar sublayer would exist. The 

changing of particle density p does not affect the curves of 

P 

Figure 28. 

A. Two-Dimensional Calculation of a Particle's 
Trajectory Under Earth Surface Conditions 

The numerical solution solves the motion of a single 
particle under various given conditions. The initial condi- 
tions calculated were: i.) particle densities of p = 3, 5 , 

P 

3 

and 7 slugs/ft , ii.) friction speed u* = 4 , 3 , 2 , and 1 

fps., and iii.) a particle diameter range from 50 - 1050 
microns in equal increments of 100 microns. This represents 



125 



Figure 28. Friction Reynolds numbers for sample ' calculation 
for Earth and Mars 
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165 separate cases calculated for each of several different 
representations of lift force as well as the case of no lift 
force with only a given initial upward velocity component. 
Functions of the lift force were developed for both the case 
of a turbulent boundary with no viscous lamxnar sublayer and 
with one that has a laminar sublayer. 

1* Case o£ 2£ lift f° rc g_ 

The dimensionless equations of motion were numerically 
solved for .an initial upward speed equal in magnitude to the 
friction velocity, i.e., 

♦ _ „ ( 6 . 1 ) 
w = z = u* 
p 

No lift force was present in these calculations and particle 
motion is due solely to the initial prescribed vertical 
velocity. For the turbulent boundary layer flow both the case 
of a laminar sublayer present and the case with no laminar 
sublayer existing were calculated. Figure 29 records the plot 
Q f t he value of maximum height (in millimeters) obtained 
during the particle's trajectory versus the particle diameter 
size D for friction speed u* = 4, 3, 2>, and 1 fps. Both 
types of turbulent boundary layer flow are presented. The 
solid line represents the case of a fully turbulent boundary 
layer while the dashed line represents the situation of a 
laminar sublayer present in the flow. This plot is for a 
constant particle density of p p = 5 {slugs/cubic ft) (the 
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average density of typical sand particles). The trends with 

respect to u* exhibited by particles of density p = 3, 7 are 

P 

similar although the maximum heights are lower for lighter 
particles and higher for the denser materials. A typical 
example is a friction speed u* = 3 fps for the fully turbulent 
boundary layer flow for a 450 micron diameter particle for 
which the maximum trajectory heights obtained by the different 
particle densities p = 3, 5, 7 were 27.8, 31.8, and 34.2 mm, 

ir 

respectively. Since the initial momentum given to the particle 
increases with increasing particle density, consequently the 
maximum height obtained in the trajectory is increased. This 
trend does not jibe with observation. The case of fully 
turbulent boundary layer flow yields a higher maximum height 
than that of a laminar sublayer flow for larger particles. 

This effect rapidly diminishes with decreasing friction speeds 
as exhibited by the coalescing of the curves for u* = 1 and 
2 fps. An apparent trend is that for a constant u* and 
increasing particle size the maximum height of the trajectory 
approaches a constant value. This trend is basically wrong 
since the trajectory’s maximum height must at some point begin 
to decrease with increasing size. At some uniquely defined 
0^ the constant friction speed will be equal to that of the 
threshold for which particles larger than this physically could 
not lift off. Since the threshold friction speed versus fric- 
tion Reynolds number curve is a double-valued curve, there is 
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also a lower limiting particle diameter below which particles 
will not lift off. The maximum difference between the two 
types of calculation (laminar sublayer, no sublayer) is never 
greater than. :4%. The result of these calculations show that 
for larger particles the initial upward (or maximum) velocity 
is not constant but should be a decreasing function of particle 
size as well as a function of friction speed. It will be 
shown later that if u* is held constant that w 0 (the maximum 
or initial vertical velocity component) is nearly constant for 
smaller size particles and then monotonically decreases for 
increasing particle size. 

2. Case of lift function for turbulent boundary layer with a 
laminar sublayer 

The analytically derived lift coefficient function (see 
Equation 4.41) from Saffman's (1965, 1968) original work was 
applied to a fully turbulent boundary layer situation and 
incorporated into the equations of motion. Particle trajec- 
tories were calculated and compared to the limited experiments] 
data available. There was no initial velocity given to the 
individual particles. The resultant trajectories were found 
to be small in comparison to experimental work performed by 
Zingg (1953). 

Zingg's data do not give the average trajectory of 
particles but rather an average height (particle mass flux 
above equals mass flux below) calculated from experimental 
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distribution of particle material caught in a vertical trayi 
An average particle size is calculated from the particle 
distribution and associated with an average height. An 
empirical equation was developed by Zingg to match the experi- 
mental data. Figure 30 shows a plot of the average height 
versus particle size. Although Zingg 1 s equation matches his 
experimental data, the trend of indefinitely increasing 
particle height with increasing particle diameter for a 
constant friction velocity is invalid. The same reasoning of 
the previous section can also be applied here and Zingg 1 s 
empirical curves predict particle motion for conditions for 
which friction speeds are far below that of threshold. This 
is obviously in error. Although Zingg’ s empirical curve 
breaks down for large particles it should be valid for the 
vicinity in which the experimental data were taken. Although 
a direct comparison of Zingg' s data to the present calculations 
is not strictly valid, the comparison is of some interest. The 
range of calculated maximum height should be from -the average 
height of Zingg' s data to approximately two times it, assuming 
the average height of the data to be somewhere between one-half 
the maximum height and the maximum height. 

The effect of the wall on particles near and resting on 
the wall was not considered in the above calculation, however, 
it appears to have an important effect on both the lift and 
drag coefficients. From Saffman (1965), on the wall effect of 
particles in simple shear flow: 



Average height , (mm) 
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The wall effect acts in two different 
ways. First, the extra drag due to the 
walls will make the particle lag behind 
the fluid; this relative velocity will 
be viscosity independent when the 
viscosity is large, and will depend only 
on the relative size of particle and the 
distance from the wall. Secondly, the 
flow field around the particle is 
altered to the presence of the walls and 
the inertial effects will differ from 
those for a particle in an unbounded 
flow, especially when the particle is 
near the walls. 

This altering of the flow around a particle will also affect 
the lift on the particle. The flow will be slower on the side 
of the particle nearest the wall and this change in the flow 
field will create a new pressure distribution on the particle. 
The pressure will be significantly higher on the under side of 
the particle thus effectively making the lift force much 
higher than normally predicted for flow in the absence of the 
wall. The lift and drag forces have the net effect of 
creating initial vertical motion, then turning the particle's 
motion to the direction of the stream. This phenomenon of a 
particle leaving the surface nearly vertically has been 
observed by the author and previously by others. It was 
assumed that the lift force for a turbulent boundary layer 
with a laminar sublayer acts only in the sublayer region. The 
laminar sublayer was then broken into two regions. The first 
region is the one nearest the wall where the wall effects are 
included in the lift force function empirically. The second 
or outer region does not include the wall effect. For Earth 
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the wall effect was empirically simulated by allowing a region 
of high lift force to act from a height of one-half the 
particle diameter to a height equal to the particle’s diameter. 
The analytic derivation of the wall effect on drag force was 
performed by Goldman et al. (1967) and O'Neill (1968). O'Neill 
(1968) obtains an exact solution of the linearized Stokes flow 
equation derived for a viscous flow about a fixed sphere in 
contact with a plane for uniform linear shear flow if the 
sphere were not present. The sphere's forces are calculated 
explicitly. Goldman et al. (1967) in an earlier paper obtained 
an exact solution to Stokes equation for the translational and 
rotational velocities of a neutrally bouyant sphere moving 
close to a plane wall for uniform linear shear flow if the 
sphere were not present. 

The empirical expressions for lift used were: 

L = 15.71pu*D^/v for i < ^ < 1 (6.2) 

P D 

and 

L = 1.615pu*D^ z/v for 1 < ^ <_ (6.3) 

' P P P * 


A velocity profile of 


u, z 

u _ * 

U* V 


(6.4) 


was used for the laminar sublayer region and a velocity 
profile of 

zu a ' 


u 1 - 
u ~ k g 


+ 5.5 


(6.5) 
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used for the turbulent portion of the boundary layer. 

Unfortunately, the laminar sublayer height lOv/u* is 
smaller than nearly all of the particle diameters calculated 
and thus the friction Reynolds number is that of flow for 
transition. In order to obtain a laminar sublayer solution 
to be later used in asymptotic matching of laminar flow to 
turbulent flow the length of the range of the second region 
was increased from lOv/D^u* to 10, i.e. 

L = 1.615 uVz/v for (1 < < 10) (6.6) 

£ - D p - 

Recalling Saffman's statement that the wall effect, depends 
only on the size of the particle and distance away from wall. 
Equation 6.6 should be approximately valid. To assume a 
constant distance that depends only on the particle diameter 
may not be correct for the fully turbulent flow where simple 
shear flow no longer exists. 

The above lift force was substituted in the equations of 
motion and again the maximum heights calculated were all found 
to be less than 1 cm high. This does not agree with Zingg's 
experimental data. In Saffman's original work the lift was in 
error by a factor of 4ir which he later corrected in 196 8. The 
above lift force in both the wall effect region and the outer 
layer of the laminar sublayer were then increased by a factor 
of 4 tt to match the lift forces predicted by Saffman's original 
work. The results of this modified lift force are given in 



135 


Figure 31. 

Figure 31 shows a plot of the calculated maximum height 

z „„ versus the particle diameter for the latter case (for 
max 

lift force in the second region of the sublayer as described 
above) . Of course these calculations are probably in error 
for the case where a laminar sublayer does not exist but will 
be utilized to approximate a solution in the transition region 
between the sublayer case and the fully turbulent one. Four 
different friction speeds, u* = 4, 3, 2, and 1 fps are shown. 
The region of small particle diameter (D^ < 300) is to be 
considered more valid than the range of larger particle 
diameter (D^ > 300). These curves also exhibit a decreasing 
z max increasing particle diameter size for a constant 

value of friction speed. They show a rapid increase in 
for a range of from 50 microns to 300 microns. 

Figure 32 shows the same type plot with the exception 
that the particle density is varied and the friction speed has 
a constant value 'of u* = 2 fps. The trend in this plot is 
for lower density particles to rise higher under identical 
conditions than those for heavier particles. This trend is 
opposite to that of the case where no lift force exists and 
force motion is created by setting w c . = u*, where the heavier 
density particles went higher under identical conditions 
(since they have more momentum) . The trend of lower density 
particles rising higher in Figure 32 seems conceptually valid. 
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Figure 33 exhibits the dimensionless maximum vertical • 

s 

velocity ratio w Q /u* versus the particle diameter for the cor- 
responding curves of Figures 31 and 32. Plausible trends are 
displayed in these plots. Here the value of initial dimension- 
less vertical velocity w e /u* decreases with increasing particle 
size which has to occur to meet threshold conditions for the 
large particles. These trends show a rather large ratio of 
w 0 /u* for larger particles (D > 400 microns) which is a 

ir 

result caused by the large extent to which lift was allowed to 
act (up to z = 10D ) , which is larger than the laminar sub^ 
layer would be if it could exist. This may, however, better 
approximate a particle initially at rest in a region of higher 
surface roughness, z 0 . 

3. Case of lift function for fully turbulent boundary layers 
The lift force and coefficient developed for a fully 
turbulent boundary layer flow in section V was substituted in 
the equations of motion and a velocity profile of 

ss ~ log(30z/D ) (6.7) 

u* k 3 p 

used, where the roughness height is 

z 0 = D /30 (6.8) 

P 

The resulting trajectories were again low in comparison 
to experimental data. The lift coefficient was increased by 
a factor 4 tt, i.e.. 


P = 5 (slugs/ft 3 ) 


200 400 600 800 

D_ ~ (microns) 


1000 


1200 


Figure 33. Maximum vertical velocity for Earth, effect of friction speed (laminar 
sublayer) 
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C T = 32.48 [(u*z/v) log (30z/D^) ] " 1/2 (6.9) 

ii p 

Figure 34 displays the resultant calculations of maximum 
height z max versus the particle size D p using the above lift 
coefficient. The calculated maximum heights are larger than 
that of the laminar sublayer flow of Figure 31. The curves 
show a sharp decrease of z max for increasing (D^ > 400) 
which must exist if threshold data are to be satisfied. These 

» 

curves (for constant u*) probably are incorrect for the smaller 
particles where a laminar sublayer would exist. The flow for 
the majority of these curves in the transition region - and the 
calculated values would not be expected to identically repre- 
sent the physical process. However, in the transition region 
especially at large friction Reynolds numbers the flow 
resembles the characteristics of turbulent flow more so than 
laminar flow. For friction Reynolds number greater than 70 
the trends the solution exhibits would be expected to be valid 
although the relative values might not be exact since neither 
analytical nor empirical values of lift coefficient in the 
transition and fully rough regions are known. 

Figure 35 shows a plot of z max versus D p for a constant 
value of u* = 2 fps for three different particle densities, 

3 

p = 3, 5, 7 slugs/ft . Again the lighter material obtains a 
P 

hiqher z than that of a heavier density particle under 
3 max 

identical conditions. This shows that lighter materials are 




Figure 34. Maximum trajectory height for Eartt 
turbulent) 



(microns) 


, effect of friction speed (fully 
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much more susceptible to gusts of wind and the lighter material 
is the first to move with an increasing wind speed. With 
decreasing densities of material under a constant set of condi- 
tions the' particle goes higher than normally would be expected 
for the denser material. 

Figure 36 records the maximum vertical components of 
velocity divided by the friction speed versus particle diameter. 
Here the trends agree with that of Figure 33 but a sharp 
decrease with increasing particle diameter is noted. This 
conceptually seems more plausible, as the larger the particle 
the slower would- it lift off under similar conditions. The 
values of the velocity ratio of larger particles (D^ > 500) 
have a range from 0.5 to 1.0 for the friction speed u* = 4. 

This is more realistic than the higher values displayed in 
Figure 34. 

4. Composite results of laminar and turbulent flows 

Although the majority of points calculated in the fully 
turbulent boundary layer are actually in the transition 
region the solutions to the particle motion are a better 
approximation than that of the laminar sublayer for reasons 
discussed earlier. With the laminar sublayer and turbulent 
flows known for the same set of initial conditions a composite 
curve of versus D may be made assuming the solution 

exists somewhere between the laminar and turbulent cases and 
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and asymptotically meets each other as flow turns from 
completely laminar to transition and from transition flow to 
fully turbulent flow. 

Figure 37 is such a plot which displays the solution for 
a fully turbulent boundary layer and one with a laminar sub- 
layer for a constant value of friction speed u* equal to 2 fps. 
Also displayed in the figure are the curves for Zingg's experi- 
mental data and for the case where no lift force acts on the 
particle, but motion is produced by giving the particle an 
initial vertical velocity component of w 0 = u*. Zingg's 
experimental points are also shown in the figure as hollow 
circle symbols. The single solid circle symbol that lies on 
the ordinate axis is the empirically determined threshold 

value and thus the curve of z„ versus particle diameter D 

max p 

must reach zero at this particle diameter. 

For lower particle diameters the laminar sublayer solution 
should be valid at least up to the particle diameter size that 
corresponds to the upper limit of the particle friction 
Reynolds number (R^) range for laminar sublayer existence. In 
this case that is a of 250 microns which has a corresponding 
Rf of 10.935 which is in the lower end of the transition range. 
At larger values of the correct solution for the transition 
region should asymptotically match the fully turbulent case at 
a R^ of approximately 70. The largest R^ obtained in this plot 
is for a of 1050 microns which corresponds to R^ = 45, thus 



(mm.) 



Dp~ (microns) 


Figure 37. Maximum trajectory heiaht for wa-v-t-v, * _ 
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the interpolated solution for the transition region would be 
expected to meet that obtained for the fully turbulent 
boundary layer at R f = 70. The dashed-dotted line is the inter- 
polated solution for the transition range that approaches the 

laminar sublayer calculation at D ^ 250 microns (correspond- 

P 

ing to Rj - 5) . It would approach the fully turbulent boundary 
layer solution at — 70 but does not reach it in this case 
since threshold conditions occur in the transition region. 

Although the final estimation of the z versus D curve 
^ max p 

lies above the experimental data points of Zingg, the solution 
appears to be valid since Zingg' s experimental points represent 
intermediate positions from one-half the maximum height to the 
maximum height as discussed earlier. Thus the estimated curve 
has the proper shape to compare favorably with experimental 
points and also to obey threshold conditions. This basic trend 

of increasing z „ for smaller D and obtaining a maximum z 

for some optimum particle size and then z max decreasing with 
increasing to the threshold condition appears to be con- 
ceptually correct. All 167 separately calculated cases have 
identical trends and same shaped curves but this particular 
case was selected to exhibit the experimental data of Zingg. 

In many of the cases calculated, unlike the above example, the 
solutions do extend into the fully turbulent region instead of 
only the transition region as above. The shape of the curve 
remains the same and continues to retain all of the above 
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characteristics. In particular, interpolated curves 'for the 
transition region. are made to asymptotically match the solution 
of the fully turbulent region. 


5. Effect of turbulence on particle trajectories 

The effect of turbulence on a particle trajectory was 
accounted for in the equations of motion and solved numerically. 
The model of turbulence used was to assume that the streamwise 
component of velocity obeyed the time averaged velocity pro- 
file of Equation 6.7 and the vertical component of velocity 
had a cyclic velocity distribution. The cyclic vertical 
fluctuating velocity had an absolute maximum value of 1.1 times 
the friction velocity, e.g.. 


W = 1.1 u* sin wt (6.10) 

where w represents the angular velocity* of the vortex eddies. 
The representation of w was taken from Panofsky and McCormick 


(1960), 


7.. 25 u* p z .+ z 0 

“ = 5 [_ 

where the roughness z Q is assumed equal to D /30. 

P 


( 6 . 11 ) 


.7.25 u. 


03 


z 


log (1 + 30z/D p ) 


( 6 . 12 ) 
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Therefore, 

"7.25 u* 

w' = 1.1 u* sin log(l + 30z/D )t (6.13) 

z p 

This expression of w‘ was incorporated in the equations of 
motion for the fully turbulent Earth boundary layer Equation 
6.9 and the motion of a particle calculated. 

The results of the calculations show that turbulence 

plays a minor role in the particles motion (D _> 100 microns). 

For example, the calculated trajectory for a particle diameter 

of 100 microns with a friction speed u* = 2 fps and p =5 

P 

3 

slug/ft shows only a 3.5% increase (with respect to the w' = 

0 calculation) in the calculated z and only a 1.1% increase 
in the final velocity of the particle u^,. Since the case of 
100 micron particles would be far more susceptible to the 
effects of turbulence than larger particles this would indi- 
cate that turbulence plays a minor role in altering the 
particle's path. An exception to this generalization is when 
the friction speed u* is of the same order of magnitude or 
greater than the terminal speed -of a particle. Turbulence 
then completely dominates the particle motion and the particle 
is then said to be in suspension. 

To extrapolate the above results to the case of Mars 
would be very hazardous since the basic structure of the 
Martian boundary layer is not known at the present time and 
also for nearly all conceivable flow conditions a rather large 
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laminar sublayer would exist. In the laminar sublayer the 
effect of turbulence is negligible and therefore the particle 
would only be affected by turbulence outside of the laminar 
sublayer. At these heights the velocity of the particle is 
relatively high in comparison to the estimated turbulent 
fluctuations and therefore need not be included in the calcu- 
lation (at the risk of an approximate 4% error). Thus, the 
effect of turbulence in the Mars calculation was not felt to 
be of importance in accurate solution of the particle path and 
was not included in the- calculations. Also the computer 
calculating time for the turbulent condition increased by a 
factor of 10 times over that for flow without simulated 
turbulence . 


B. Two-dimensional Calculation of Particle 
Trajectories under Mars Surface Conditions' 

The numerical solutions obtained for the Mars surface 
conditions were performed in a similar fashion as those cal- 
culated for Earth with the exception of a fully turbulent 
boundary layer case which was not calculated. All cases 
calculated for Mars have a friction Reynolds number less 
than seven and a laminar sublayer exists for all of these 
cases. It was necessary to consider slip flow around the 
particle caused by the rarefied atmosphere of Mars. The 
numerical coefficients used for these calculations as well as 
the different variables used in the calculations are given in 
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Table 3. 

The values of friction speed u* were selected by approxi- 
mately matching the ratio of flow friction speed to that of 
threshold, u*/u* t , to that of the cases calculated for Earth. 
These calculated values of u* t were' taken from the work of 
Iversen et al. (see Section III) and are presented in Table 4. 
It was assumed that a similar range of particle densities 
existed on Mars as on Earth, as well as a similar range of 
particle size distribution. The atmospheric density was taken 
from the work Iversen et al. (1975a) . 

1. Case of no lift force 

The dimensionless equations of motion were solved for an 

initial upward vertical component of velocity equal to one- 

tenth and one-twentieth of the prevailing friction speed u*. 

These values of initial upward velocity were calculated from 

trajectory calculations with lift and from the work of Bagnold 

(1941) in which he assumes the quantity (u p w 0 /gy ) is a 

r max 

constant for any universal set of conditions. These values of 
w 0 were also supported by numerical calculations for the 
laminar sublayer case. 

Figure 38 displays a plot of z max versus D^, for the no 
lift force case with w 0 set equal to one-tenth u*. Similar 
trends as those calculated for Earth with, no lift force present 
are observed. Here the turbulent boundary layer with a laminar 
sublayer velocity profile was used. This plot is for a 
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Table 4. Threshold friction speeds (fps) 




Earth 



Mars 


p p (slug/ft 3 ) 3 

5 

7 

3 

5 

7 

D (microns) 
P 






50 

0.531 

0.686 

0.700 

9.024 

9.578 

10.095 

150 

0.540 

0.698 

0.826 

7.464 

8.757 

9.808 

250 

0.698 

0.901 

1.066 

7.666 

9.164 

10.344 

350 

0.826 

1.066 

1.261 

• 7.93 

9.60 

10.865 

450 

0.9 36 

1.209 

1.430 

8.23 

10.05 

11.445 

550 

1.035 

1.336 

1.581 

8.53 

10.465 

12.02 

650 

1.125 

1.453 

1.719 

8.78 

10.85 

12.56 

750 

1.209 

1.560 

1.846 

9.04 

11.22 

13.11 

850 

1.287 

1.661 

. 1.965 

9.29 

11.60 

13.6-1 

950 

1.36 

1.756 

2.078 

9.52 

11.93 

14.06 

1050 

1.43 

1.846 

2.184 

9.75 

12.24 

14.48 

constant 

. particle density 

of p p ■ 5 

(slugs/ft 

: 3 >. 


It 

might be 

noted that similar 

trends to 

those 

discussed 

earlier 

(Earth case of no 

lift, w 0 = 

= u* ) for 

changing particle 

density 

also are 

present in these calculations and the same 

analysis 

can be 

applied to 

' these as 

was done 

to the 

Earth case 


Recall that the main point made earlier was that the trend of 





max 



Figure 38. Maximum trajectory height for Mars, no lift, w 0 = 0.1 u 
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heavier density particles lifting off higher than those of 
lighter density for identical conditions is incorrect. The 
Wo = constant calculations are thus used only to gain insight 
to the overall transport process. 1 1‘ might be noted that in 
comparison to Figure 29 (the Earth case) that the Earth curves 
fall off much more quickly at the smaller particle size than 
that of Mars. An apparent explanation of this phenomenon is 
the slip flow effect of Mars. Even though the w 0 on Mars is 
much smaller (on the order of 10% of Earth's) the shape of the 
calculated curves and their trends should be similar to 
Earth with the exception of the effect of slip flow. 

The effect of slip flow can be observed in Figure 30, 
where the plot shows Z max / D p versus for p p equal to 5 
(slugs/ft^) and w 0 = 0.1 u* (with no lift force) for a 
friction velocity u* = 25 fps. The solid line re pres ents the 
solution of the equations of motion including the slip flow 
effect of C. N. Davies (1945). The dashed line displays the 
solution for continuum flow (no rarefied gas effects included) . 
There appears to be very little difference in the solution 
down to a particle diameter of approximately 250 microns. The 
effect of slip flow drastically increases for 50 micron 
diameter particles. This effect is 'identical for all cases of 
^ars calculated since the effect of slip flow is independent 
of the lift mechanism employed to particle motion. Therefore, 
it will suffice to examine closely only one case of the effect 




D ~ (microns) 

P 

Figure 39. Maximum trajectory height for Mars, effect of slip flow in a rarefied 
atmosphere 
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of slip flow compared to continuum flow as was done here. 

Figure 40 shows the relation between Z max / D p versus 
curves for w 0 = 0.1 u* and w 0 = 0.05 u*. Each trajectory was 
calculated with identical conditions (with the exception of 
initial velocity) which were u* = 25 fps and p = 5.00 slugs/ 

ir 

ft^. The figure displays similar trends and only the magnitude 
of the curve is increased for increasing initial velocity. As 
the initial momentum is increased, the particle maximum energy 
or height will also increase. The validity of one curve over 
the other curve in Figure 40 remains unknown but a numerical 
investigation of the case of a lift present in the laminar 
sub.layer region with the initial vertical and horizontal 
velocity set equal to zero will indicate which curve in Figure 
40 more closely represents the trends exhibited by the solution 
of the equations of motion for a lift force present.' 

2. Case of lift function for turbulent boundary layer with a 
laminar sublayer 

The analytically derived lift coefficient function (see 

Equation 4.41) based on Saffman's (1965, 1968) work was used 

in the solution of the equations of motion. The solutions (as 

in the Earth case) of the particle's motion resulted in very 

low z and are not believed by the author to be valid, 
max 

Consequently, the factor of 4r was applied to the lift 
coefficient and the equations solved again. This is the same 
process or technique that was used in the Earth calculations 




Figure 40. Maximum trajectory height for Mars, no’ lift, w 0 


0.1 u* and w Q 


0.05 u 
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and is believed to be as valid for these calculations (Mars) 
as it was for the Earth. Essentially this states that the 
solutions obtained in the Mars calculation should be of the 
same order of accuracy as the solutions of the Earth case. 
Recall that the Earth solution was definitely within an order 
of magnitude of related experimental data and believed to have 
a valid and accurate trend in analysis of the experimental 
data of Zingg (1953) and supporting threshold work of Iversen 
et al. (1975a) . Therefore, extending this analysis, the Mars 
calculations should be of the same order of magnitude as the 
physical situation and quite possibly would describe the 
particle motion very accurately. Of course, these claims can 
not yet be supported by experimental data or visual observa- 
tion. 

A velocity profile for the viscous laminar sublayer and 
the turbulent portion of' the boundary layer were the same as 
Equations 6.4 and 6.5, with the appropriate values of Mars 
used for the various constants. 

Again as for Earth, the sublayer was broken into two 
regions; one accounting for the wall effect where the lift 
coefficient was 


C T = 41.12 D v/u*z 2 for i < z/D < 1 
L p * Z — p — 


(6.14) 


and for the second region was 


1 < z/D < 


lOv 


“ P ~ u *° r 


C_ = 12.91V/TT AV D 

Jj jO p 


for 


(6.15) 
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where both of the above are multiplied by 4 tt. The fully 

developed equations of motion are presented in Appendix B. 

Figure 41 shows the solution of the equations of motion 

using the above assumptions. This figure displays z max versus 

3 

D for p equal to 5 (slugs/ft ) and u* = 25, 20, 15, and 10 

Ir ir 

fps. The z max calculated here are generally higher than those 

calculated for Earth recalling that the ratio of u*/u* t is 

approximately equal thus making a comparison valid. The 

trends of the Mars, calculations are basically different from 

those of the Earth calculations and this altering of the 

curves (for D > 250 microns) are caused by boundary layer 
P 

profile differences due to different particle friction Reynolds 

number R^. However, no physical explanation of the shape of 

the curves for larger particles (D^ > 250 microns) is offered. 

The point of interest is the gradually decreasing slope of the 

curves for larger particles for decreasing friction speed u*. 

It is apparent that these curves (for u* — 15 fps) must be 

returned to the ordinate axis with a further increase in 

in order to meet threshold conditions predicted by Iversen 

et al. (1975a) . Thus there should be a double peak curve of 

z , versus D . The predicted threshold for 450 microns is 
max, p c 

10 fps so the lower curve should appear closer to the dashed 
line as shown. 

An alternate way of presenting the curves of Figure 41 is 
to normalize the maximum height z by the particle diameter 



10 


Figure 41. 



Maximum trajectory height for Mars, 
sublayer) 


effect of friction speed (laminar 
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D p' e ‘ g * Z max /D p versus D p - Figure 42 is such a plot of the 
identical data as presented in Figure 41. The unusual char- . 
acteristics of the curves in Figure 41 disappear and a unique 
set of curves results. Although the set of curves presented 
xn Figure 41 are nicely defined they may be misleading since 
the normalizing factor D p is continually changing which makes 
relative comparison difficult between different size particle 
diameters D p . The dashed line shown corresponds to that of 
Figure 4-1. 

Figure 43 records a plot of z max versus D p for a constant 

u* of 20 fps and three different particle densities o = 3 5 

f r 

and 7 slugs/ft 3 . The trends of these plots are opposite to 

those of the case of no lift force. The curves of Figure 43 

infer that light density material lift off to a higher z 

' 3 max 

than heavier density material under the same flow conditions. 
These trends agree with those of the Earth case calculation 
with a laminar sublayer. The relative increase of z max with 

decrease of is similar to the trends exhibited in the Earth 
calculations . 

Figure 44 shows a plot of w 0 /u A for the corresponding 
conditions of Figure 42 and 41. It shows a nearly constant 

w °/ u * 0*1 for all friction speeds u* for D greater than 

’ ]? 

450 microns. An interesting point on the curve is' at D = 

P 

550 mxcrons where all w Q /u* ratios are almost identical. Here 
all the curves reverse position in relation with each other. 
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Figure 43. Maximum trajectory height for Mars, effect of particle density (laminar 
sublayer) 
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Figure 44. Maximum vertical velocity for Mars, effect of 
friction speed (laminar sublayer) 
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At smaller (D < 350 microns) the ratio w 0 /u* increases 

rapidly as is the case of the fully turbulent boundary layer 

calculation for Earth. .This exhibits the susceptibility of 

the smaller particles to movement by wind. 

Figure 45 shows a plot of z_ a „ versus u* for four dif- 

max x 

ferent particle diameters D = 250, 550, 750, and 1050 microns. 

P 

An interesting point of this plot is that the variation of 

z max ^ or a cons 't : ant u* for a large range in D is relatively 

P 

small and all seem to cross for a u^ range of 16-21 fps and 
reverse relative positions between the various particle 
diameters. The dashed lines connect the predicted threshold 
value (Iversen, 1975b) to curves of the numerical results. 

Several additional mechanisms of lift function were also 
calculated. Of these, four cases were calculated in which the 
effect of the wall was increased from one particle diameter to 
2, 3, 4, and 5 particle diameters. The resultant solution of 
the equations of motion shows verv little increase in the z 

max 

calculation. Both the Saffman lift and the modified Saffman 
lift (4ir times the Saffman lift) were calculated with similar 
results . 

One additional case calculated was that of the one 
earlier presented for Figure 41-45 with the exception that the 
lift function of 4ir times Equation 6.6 was allowed to act for 
a distance of 10 particle diameters. This essentially 
curtailed the region of active lift force for smaller D 

P 



max 



Figure 45. Maximum trajectory height for Mars, effect of 
particle size (laminar sublayer) 
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(D < 200 microns) and enhanced the net lift force of larger 
P 

D (D > 200 microns) . This* did not change the lift force 
P P 

produced in the wall effect region but only altered the -length 

of the outer region where no wall effects were considered. 

This type of calculation has the same type of validity as the 

similar case presented for the Earth case (see the Earth case 

of a laminar sublayer for a detailed discussion) . Figure 46 

displays a comparison of z max versus for u* = 20, and 25 

fps for the lift function where the outer lift region was 

allowed to act to 10 particle diameters and to the sublayer 

height for a constant - 7 (slugs/ft^). As can be observed 

from the figure the effect of allowing lift to act to 10 

particle diameters for the larger particles enhances the z max 

by a factor of 2 (D p > 250) and suppresses the z max for smaller 

particle diameter (D < 150) . The validity of the sublayer 

calculation is believed by the author to be more accurate than 

the comparative lift model of 10 particle diameters; however, 

the true solution can only be guessed at this time. Figure 47 

shows the corresponding curves of Figure 46 for w Q /u* versus 

D , and again similar trends to those of Figure 46 manifest 
P 

themselves . 

3. Particle rebounding at Martian surface 

An important process of particle saltation motion is the 
bouncing interaction that occurs at the end of a particle 
trajectory. This process has long been observed and it is 
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Figure 46. Maximum trajectory height for Mars, cc 




Figure 47. Maximum vertical velocity for Mars, 


D case u *~ 

d — U-Q o a o 

□ □ □ a □ 



Sublayer case 


i i 

800 1000 

crons) 


of lift model 



169 


known that the collision of a single particle with a surface 
of random composite of similar particles results in a highly 
inelastic collision. Generally the single particle in motion 
imparts the majority of its momentum to the surrounding 
particles in contact with it at the surface which, according 
to Allen (1968) , initiates several other particles in motion. 
The ensuing initial motion is nearly always vertical. These 
particles then turn and travel in low trajectories at veloci- 
ties nearly parallel to the surface at the end of the trajec 
tory.. When particles strike the surface, a particle may either 
rebound off the surface from a stationary particle or group of 
particles, or impact into the surface of particles when several 
other particles' motion may be initiated as discussed earlier. 
In the case of the rebounding particle (or for that matter 
particles initiated with a given percentage of the original 
particle's momentum) the lifting process is a combination of 
the initial upward component of vertical as well as the lift 
force function . 

Figure 48 is a plot of the normalized maximum height of 

the trajectory and the final (or maximum) downstream distance, 

z /D and y /D respectively, versus the percentage of 
max 7 p J max p 

momentum retained by the rebounding particle or which may be 
considered as the percentage of momentum imparted to a similar 
size particle upon collision. The lift function used in the 
calculations is that of 4ir times Equations 6. 2. and 6.4 for a 
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Figure 48. Effect of rebound momentum of trajectory height 
and length 
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100 and 500 microns particle diameter, however, for the case 
of the 100 microns particle the range of the lift function is 
the same as that of Equation 6.3. Nevertheless ," the normalized 
quantities should be the same for similar type lifting func- 
tions as in the case here. As can be observed from these plots 

one sees that z _ /D increases much more rapidly than Y^.,,./ 0 ™ 
max p . max. p 

for both the cases of 100 microns and 500 microns. 

There is a rapid increase in z /D for small changes in 

III ci A p 

the percentage of momentum. An unusual characteristic of the 

100 microns y /D curve in the figure where a 4-5% momentum 
max p 

of original impact the Y max /Dp takes an unpredicted decrease 

in value. The value of curves is unity at zero momentum 

transfer and is an ordinary lift function calculation. 

In connection with the rebounding particle it might be 

of interest to note that the particle returning to the surface 

seldom makes an angle with the surface of greater than 3 

degrees and is often less than one degree. Th-is value is 

unusually low but is characteristic of all Martian trajectory 

calculations . The typical collision for Earth conditions , as 

reported by Allen (1968), is in a range of 5 to 10 degrees 

which agrees with the Earth calculations. An important fact 

noted here then is that the Martian trajectories are much more 

parallel to the surface at collision than Earth trajectories. 

As a result the erosional effect may be greater on Mars than 

on Earth. The y m , v /z m=v ratio is much greater on. Mars in 

max max 
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comparison to Earth calculations. In a paper by Andres (1969) 
the Mars particles trajectories were calculated and the surface 
interactions were taken into account. From Andres: "It was 

assumed that the incident particle upon impact with the surface 
loses the velocity component normal to the surface; that the 
particle rotation and remaining tangential velocity component 
adjust to the condition of rolling contact; and that the 
particle leaves the surface with these adjusted tangential and 
angular velocities in the direction given by the local surface 
slope." For the cases presented in Figure 48 this would infer 
that less than 1% of the impact momentum would be lost since 
the collision angle is very small. This assumption of Andres 
seems to be basically wrong and the actual situation appears 
to be a highly complicated interaction of several particles, 
none of which receive a large percentage of momentum - . 

Other numerical calculations of particle trajectories 
under Mars surface conditions were performed by Henry (1975) 
in which particle impact on wind sensors was investigated. 

C. Simulated Flow Around a Crater 

An interesting application of particulate flow is that of- 
simulating particle motion in the downstream wake of crater 
flow-. A detailed discussion of the characteristics of the flow 
field around a crater are given in Section V and pictorially 
displayed in Figure 3. The main aspect on the downstream wake 
flow of a crater is the existence of two counterrotating 
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vorticies formed from the horseshoe vortex which is wrapped 
around the front and sides of the outside crater rim. The 
direction of .rotation of the two downstream vorticies is the 
same as the trailing vorticies created by a finite airfoil 
moving through a fluid. 

A set of three- dimens xonax equations or motion were 
developed and' programmed with capabilities of analytically 
input crossflow functions for the v and w components of 
velocity as well as the stream component u (recall the 
c.oordinate used here is the same as that of Figure 9). To 
numerically simulate the downstream flow of a crater with a 
vortex present the flow is assumed to be symmetric and stable 
around the crater. Thus it is only necessary to solve one- 
half the flow field (separated by the axial center x - 0). 

A Rankine vortex is used in the numerical calculations 
to form the crossflow. From visual observation and studying 
the details of Tani's et al. (1962 > , Gregory and Walker's 
(1951), and Sedney's (1973) paper, the single vortex core is 
positioned at a height z 0 = 2 cm above the surface and a 
transverse distance of x„ = 6 in. ( £ 1-5 cm) off of the axial 
centerline. The crossflow component of velocity V Q is then 
expressed as 

V. = — - r - 9 for r <_ r 

2**1 


(6.16a) 
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and 

V Q = r/2Trr for r > r i (6.16b) 

.where Y is the strength of the circulation, r^ is the radius 
of the vortex core, and r is the distance of the particle from 
the center of the vortex core, i.e., 

r - / (z-z 0 ) 2 + (x-x Q ) 2 (6.17) 

The value of the strength of the vortex T was selected 
from the experimental results of Section V, and given a value 
that corresponds to a maximum crossflow speed (at r = r^) of 
10% of freestream value, i.e., 

V Q = 0.1 u m (at r = r ± ) (6.18) 

And from experimental results it is found that 

u ot = 23 u* (6.19) 

and thus substituting this into the above equation yields 

V 0 = 2.3 u* ' (at r = r^ (6.20) 

or 

r = 14.45 u* r ± (6.21) 

and for r^ = 1 cm this is 

T = 0.47396 u*[ft 2 /sec] 


( 6 . 22 ) 
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This expression predicts a velocity flow through the 
surface which physically can not exist and therefore was 
corrected to account for this wall effect. The empirically 
approximated wall effect is assumed; to have the same funda- 
mental properties as the function of axial velocity. This 
states that the wall effect of the vortex motion tends to 
curtail the flow in much a similar matter as the axial velocity 
profile, i.e. 

V fi = f, (u* r, ) f, (u/u 0 ) (6.23) 

eff 1 ' 1 2 

j 

where V Q is the effective crossflow velocity, f, is the 
0 eff ' 1 

unbound vortex motion function, and f 2 is the empirical 

function of wall effect. A suitable f 2 was selected as 

f 2 = (u/u 0 )°* 15 (6.24) 

where u e is a reference velocity calculated at the height of 
the center of the vortex core (2 cm). Thus, at a height of 
the vortex core the motion of the vortex is unaffected by the 
wall. 

The components of crossflow v and w then are expressed as 

v » V (z 0 -z)/r (6.25) 

eff 

w = y (x-x 0 )/r (6.26) 

• eff 

for the selected direction of the vortex. 
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The equations of motion including the crossflow given 

above were solved for the Martian surface condition of 

3 

p =5 (slugs/ft ) and a friction velocity of u* = 20 fps 
P 

with a lifting function of Equations 6.2 and 6.3 for both a 
100 and 500 micron particle diameter. Generally speaking the 
motion of saltating particles is severely inhibited on the 
downwind side of the vortex motion and greatly enhanced on 
the upwind side. This causes increases in the final velocity 
of the particle's trajectories, u_, and also in y„,^ , z mri , 
and the maximum upward component of vertical velocity, w„. 

Figures 49 and 50 show plots of the above quantities 
(normalized by the vortex-free solution) versus the transverse 
distance x away from the axial centerline for both a 100 and 
500 micron particle diameter case, respectively. The location 
of the vortex position is denoted by a single vertical dashed 
line. 

Figure 49 has two different vertical scales, the left one 
for w and V_ and the right one for the and y„ . The 

plot displays the enormous effect the vortex crossflow has on 
the particle's trajectories. For transverse length less than 
5 inches the vortex's vertical component velocity is downward 
and inhibits the particle motion in all four variables 
displayed, but for larger values of x (x :> 5) the motion is 
tremendously increased in all variables. Both x max and z max 
obtain values nearly 50 times their normal value at the vortex 
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Figure 49. 


Transverse length x , (in) 

Effect of initial particle location with respect to vortex on traiectorv 
characteristics (D p = 100u) J - y 




Transverse length x , (in) 


Figure 50. Effect of initial particle location with respect 
to vortex on trajectory characteristics 
D — 
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location and then gradually return to unity with increasing 

distance. The normalized final velocity u^, and normalized 

maximum vertical velocity w Q obtain values of 3.63 and 8.8, 

respectively. An interesting point is that the Up does not 

immediately return to unity but retains a value nearly 

constant of 3.5. This would decrease the collision angle of 

the surface compared to normal trajectories without vortex 

presence. Also as the transverse length is increased (x >_ 5 

inches) the z decreases but the u„ stays the same. This 
max ,r 

would also decrease the collision angle further with increase 
in transverse length. The net effect of this would be to 
increase the erosion on the outward side of the vortex 
transverse location. These numerical results seem to support 
the theories of a similar nature expressed by and the experi- 
mental results of Greeley et al. (197 4) and Iversen 'et al. 
(1975a), This is a remarkable result which also supports 
imagery received from Mariner 9 _ (see Figures 1 and 2) . 

Figure 50 shows the similar case for a 500 micron 
particle diameter but the effects of the vortex motion are 
considerably less. 

Figure 51 displays the lateral motion (or spiral) of a 
single particle initially started at several various transverse 
positions. All of the trajectories show a motion of the 
particle that generally follows the direction of the vortex 
velocity. These trajectories are for the corresponding 
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Figure 51. Lateral and vertical particle position with respect to the vortex. 
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calculations of Figure 49. These curves show that the smaller 
particles are greatly affected by the vortex motion and obtain 
heights of two orders of magnitude higher than normal. It 
also infers from the curvature of the particle trajectories 
that there would be considerably more particle interaction 
(collisions) in the vicinity of the vortex cores. 

The net conclusion of these vortex calculations is that 
the final velocity of the particle u^, is increased on the 
outer edge of the vortex core thus causing a greater increase 
in surface erosion rate. Also trajectory calculations show a 
curvature of the particle motion that follows the vortex motion 
which affects smaller particles much more than. larger 
particles . 


C > ^ 
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VII. CONCLUDING REMARKS 

The experimental pitot tube study of the three-dimensional 
wake flow field around the model crater performed in an 
environmental wind tunnel proved to be very important. A 
qualitative description of the surface shear stress was ob- 
tained which exhibits large shear stresses immediately down- 
stream of the crater centerline. The data also displayed 
trends that occurred in the dimensionless velocity profile 
with different locations in the flow' field downstream of -the 
crater including velocity defects possibly resulting from the 
horseshoe vortex system. In addition to the pitot probe, a 
three-dimensional hot- film sensor was used' to measure magni- 
tudes of crossflow. These were used later in three- 
dimensional numerical crater simulation flow. 

For the two-dimensional numerical study of particle flow 
the importance of the lift force in the equation of motion was 
demonstrated. Empirical lift functions were developed for 
both laminar sublayer and fully turbulent flow for Earth (that 
approximates the limited experimental data) and was used to 
approximate transition flow. The effect of an idealized 
cyclic numerical turbulence model showed turbulence to be a 
minor factor in the determination of particle trajectories 
(D > 100 microns) . 

IT 
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The empirical lift functions were used to calculate 
particle trajectories for the Mars surface conditions and 
should be as accurate as for Earth. The effect of slip flow 
was shown to be an important factor in small diameter particle 
trajectories (D^ < 250 microns) and therefore was included in 
^11 Mars calculations. In comparison to Earth trajectories 
for the same friction velocity ratio, u*/u* t , the Mars 
trajectories exhibit higher z max ' s (up to approximately 40%) 
and longer y 's than that of Earth. The ratio was 

considerably larger on Mars and exhibited a much lower surface 
collision angle (generally less than 3°) than Earth's (usually 
from 5° to 10°). Some other significant results are that the 
maximum upward velocity ratio w 0 /u* is approximately one-tenth 
of that of Earth's, and the final velocity u^, is much hiqher 
on Mars thus causing a high erosional effect. 

The effect of momentum exchange of inelastic collisions 
of particles at the surface for Mars was investigated and it 
is found that smaller particles rebounded several times higher 
than normally would be predicted, while retaining after col- 
lision only a small percentage of its total momentum. This 
effect decreased with increasing particle size. 

Finally, a combination of the experimental data and two- 
dimensional particle calculation was used in the three- 
dimensional particle flow in the wake of a crater. The 
significant results of this combination of a turbulent layer 
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and a vortex flow include curtailment of particle motion on 
the side of the vortex closest to the wake centerline and an 
enhancement of the trajectories on the outer -side of the 
vortex. An increased .value of u p and a decreased value of the 
collision angle compared to vortex free motion result both of 
which add to the erosional effect.;' Also the small particle 
(0^ < 200 microns) trajectories were highly affected by the 
vortex xuotion. 

However, the results of this study axe by no means con- 
clusive and would indicate the need for future work* Of the 
future work an interesting research area would be the experi- 
. mental investigation of the actual lift force exerted on small 
particles and also of interest would be an experimental 
investigation of the average paths of particle motion under 
various wind tunnel conditions. This could, possibly be 
accomplished by the use of high speed movies. 

Numerically, the ultimate goal of this study would be to 
combine a three-dimensional flow field with the particle equa- 
tions of motion to study the motion of a particle in flow 
around a crater. Several particles could simultaneously be 
followed through the flow field. The specific application of 
this process to Martian conditions could possibly match the 
imagery (received from Mariner 9j of what is believed to be 
caused by flow over Martian craters provided the movement of 
many particles is allowed to occur. The resultant type of 
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program could also be used to test .newly designed structures 
that might be subject to snow, dust and sand storms. This 
ultimate goal is quite complex and would involve a large amount 
of computer storage and computer time to solve for the 
trajectories of hundreds and even thousands of particles in 
simultaneous motion. Eventually, perhaps a dual process could 
be involved where the particle equations of motion are solved 
numerically and then the partial differential equations of 
motion could be solved thus yielding a new flow field. This 
process could be repeated until the problem is solved. Cur- 
rently, most such solutions are beyond the storage capacity of 
present day computers. 
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X. APPENDIX A 

Listed here are the formulas used to calculate the drag 
coefficient numerically in the computer programs that numerical 
solved the system of ordinary differential equations of the 
motion of the particle's trajectory. 

= 24.0/R for R £ 0.1 

C D = 22.73/R + 0.0903/R 2 + 3.69 for 0.1 < R £ 1.0 

C D = 29 . 1667/R - 3.8889/R 2 + 1.222 for 1.0 < R £ 10.0 

C D = -46. 5/R - 116.67/R 2 + 0.6167 for 10.0 < R < 100.0 

C Q = 98.33/R - 2778. 0/R 2 + 0.3644 for 100.0 < R £ 1000.0 

C D = 148.62/R - 4.75 x 10 4 /R 2 + 0.357 for 1000.0 < R < 5000.0 

C D =-490 . 546/R + 57.87 x 10 4 /R 2 + 0.46 

for 5000.0 < R £ 10000.0 
C D =-1662. 5/R + 5.4167 x 10 6 /R 2 + 0.5191 

for 10000.0 < R £ 50000.0 

C D = 0.4 for R > 50000.0 
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XI. APPENDIX B 

Two-dimensional Equations of Motion 
m p x = J pAV r Dp [ (u-x)C D + (z - w)C L ] 

m p z = 2- pAV r Dp[ (w-z)C D + (n-x)C L I-m p g 

where AV r = /(u-x) 2 + (w-z) 2 


m 


P 


=ffpD 3 /6 

P 


Dimensionless Equations of Motion 


x = 0.75 A V [ (z-w) C T + (u-x)C n l 

D I; Ju u 


i = 0.75,2- . AV [ (u-x)C_ + (w-z)C ]-g 

p p r L - ° 


where 


A A 

= x/D ; x = x/u* ; x = xD /u* 


g = gD /u* ; u = u/u* ; AV = AV /n 4 


Three-dimensional Equations of Motion 
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i pAV r D p 
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C D + L x C L 


it AT7 2_2 
V = 8 pAV r D p 


^ S + Vl 
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(B. 2) 
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(B. 4 ) 
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(B . 6 ) 

(B . 7) 

(B. 8) 
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where 
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y AL 
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z AL 


AL = (w-z) [v(w-z) -w (v-y) ]- (u-x) [u(v-y)-v(u-x) ] } 

+ { (v-y) [u(v-y) -v(u-x) ]- (w-z) [w (u-x) -u(w-z) ] } 


z — 
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